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ABSTRACT 


Satellite  x-ray  observations  indicate  that  most  solar  flares 


produce  a hot,  ^ 10  K,  dense  plasma  high  up  in  the  sun's  corona 


We  believe  that  this  gas  originates  from  lower  regions  of  the  atmos 


phere,  specifically  the  chromosphere,  where  it  is  heated  by  a flare 


arid  evaporated  into  the  corona.  Optical  observations  have  shown  that 


the  post-flare  plasma  does  not  cool  uniformly  rather,  small  cold 


condensations  form  at  the  top  of  magnetic  field  arches  while  the  bulk 


of  the  plasma  remains  at  high  temperature.  Loop  Prominence  Systems 


are  an  example  of  this  phenomenon.  Dynamical  processes  connected  with 


the  flare  event  have  been  proposed  as  a mechanism  for  these  types  of 


active  prominences 


We  have  investigated  the  cooling  of  post-flare  plasmas,  and 


attempted  to  explain  the  formation  of  loop  prominences  as  due  to  a 


thermal  instability.  At  temperatures  between  lO'^  to  10  K,  the 


solar  plasma  is  known  to  be  unstable  to  thermal  perturbations  because 


the  radiation  losses  of  this  gas  increase  with  decreasing  temperature 


A simple  perturbation  analysis  is  insufficient  to  account  for  the 


evolution  of  the  solar  plasma,  because  observed  differences  of  temper 


ature  and  density  exceed  an  order  of  magnitude.  Additionally,  the 


plasma  is  always  in  a non-equilibrium  state,  so  that,  we  must  solve 


numerically  the  equations  of  motion  and  heat  transfer 


We  have  developed  a one-dimensional  model  for  active  loop 


prominences.  Magnetic  fields  present  in  post-flare  regions  are  strong 


enough  to  dominate  the  plasma;  hence,  only  motion  and  heat  fluxes 


parallel  to  the  field  need  to  be  considered.  The  relevant  size  scales 
and  time  scales  are  such  that  single-fluid  MHD  equations  are  valid. 

We  have  Included  in  the  model  the  effects  of  gravity,  the  geomt l ry  of 
the  field  and  conduction  losses  to  the  chromosphere.  A computer  code 
for  the  solution  of  our  set  of  equations  has  been  constructed.  Basically, 
we  treat  the  system  as  an  Initial  value  problem  (with  certain  boundary 
conditions  at  the  chromosphere-corona  transition  region),  and  use  a 
two-step  time  differencing  scheme. 

Our  calculations  Indicate  that  the  thermal  Instability  mechanism 
is,  by  Itself,  sufficient  to  account  for  the  behavior  of  active  loop 
prominences.  Under  certain  conditions.  Initial  perturbations  In  the 
temperature  and  density  profiles  of  small  amplitude  (&  ^io)  and  large 
size  scale  2 X 10^  cm.)  can  grow  Into  condensations  with  temperature 

and  density  differences  of  over  an  order  of  magnitude  and  size  scales 

g 

of  less  than  10  cm.  In  agreement  with  observations,  the  conditions 
that  must  be  satisfied  are  such  that  Loop  Prominence  Systems  are  likely 
to  occur  only  In  large  flares.  The  velocities,  densities,  and  lifetimes 
that  we  obtain  for  the  loop  material  are  also  In  agreement  with 
observations . 

From  our  results  we  conclude  that  the  non-uniform  cooling  of  the 
post-flare  corona  can  be  understood  as  a direct  consequence  of  the 
temperature  and  density  dependence  of  the  radiative  losses  from  a 
high- temperature  solar  plasma. 
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1 . THERMAL  INSTABILITY 

1 . 1 Introduction 

Parker  (1953)  first  introduced  the  concept  of  thermal  instability 

if. 

to  explain  the  formation  of  cold  condensations,  T 10  K,  in  the  hot 
solar  corona,  T 10^  K (Tandberg-Hanssen,  1967)-  He  examined  the 
stability  of  a hydrogen  plasma  to  a thermal  perturbation  and  deduced 
that  if  the  total  energy  losses  minus  energy  gains  decrease  with  T, 
then  the  olasma  would  be  unstable.  This  is  because  a region  slightly 
cooler  than  its  surroundings  will  lose  energy  faster  than  the  surround- 
ings and,  therefore,  temperature  differences  will  increase. 

Parker  also  noted  that  there  is  a temperature  range  over  which 
the  radiation  losses  of  an  optically  thin  hydrogen  plasma  decrease 

with  temperature.  This  is  the  temperature  range  over  which  the  state 

4 5 

of  ionization  of  hydrogen  changes  rapidly,  10  < T < 10  . The  emission 

is  dominated  by  line  and  recombination  radiation  for  these  temperatures. 
Above  this  range,  the  emission  is  due  mainly  to  free-free  transitions; 
hence. 


(1.1) 


where  f is  the  radiative  loss  rate  (Karzas  and  Latter,  I96I).  Below 

k 

10  K,  the  mean  kinetic  energy  per  particle  is  insufficient  to  excite 
neutral  hydrogen  above  its  ground  state,  and  the  emission  decreases 
exponentially  with  decreasing  T. 

Defouw  (1970b)  has  calculated  the  radiation  from  a model  hydrogen 


plasma.  He  concludes  that  hydrogen  is  thermally  unstable  above  17,500  K. 
Thomas  and  Athay  (I96I)  find  that  the  solar  chromosphere  is  unstable  above 
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12,000  K.  Cox  and  Tucker  (I969)  (Figure  1.1),  Cox  and  Daltabuit  (I97I1, 
Tucker  and  Koren  (197I)  and  Raymond  et  al.  (I976)  have  calculated  the 
total  radiation  losses  (including  heavy  elements),  for  the  solar  coronal 
plasma.  From  Figure  1.1  we  note  that  there  is  a region,  10^  T ^ 10^, 
over  which  J*  decreases  with  T;  hence,  thermal  instability  is  possible. 

The  original  treatment  of  thermal  instability  by  Parker  (1953) 
was  faulty  because  he  did  not  take  into  account  the  dynamics  of  the 
plasma.  Field  (I965)  in  his  excellent  paper  on  the  subject  thoroughly 
analyzed  the  linear  growth  of  thermal  perturbations  in  a fully-ionized 
plasma.  He  not  only  included  the  dynamics,  but  also  analyzed  in  detail 
effects  due  to  conduction,  magnetic  fields,  gravity,  expansion  and  rota- 
tion. Goldsmith  (1970)  and  Defouw  (1970a, b,c)  extended  Field's  results 
to  include  ionization  effects  and  Hunter  (1970)  investigated  thermal 
instability  in  a plasma  with  an  arbitrary  velocity  profile. 

We  will  deal  only  with  thermal  instability  in  post-flare  coronal 
loops.  The  important  effects  in  this  case  are  due  to  conduction  and 
magnetic  fields.  Since  the  field  dominates  the  plasma,  we  need  only 
consider  the  growth  of  instabilities  in  one  dimension,  along  the  field 
(Section  3«^)»  Field  has  derived  the  growth  rates  for  plasma  in  a 
uniform  magnetic  field;  however,  solar  fields  are  far  from  uniform 
(Harvey  et  al.,  1972;  Rust  and  Bar,  1973).  This  means  that  a constant 
flux  tube  in  the  solar  corona  will  have  a large  change  in  cross- 
sectional  area  along  its  length.  Hence,  we  derive  below  the  growth 
equation  (dispersion  relation)  for  thermal  perturbation  in  a loop  of 
varying  area,  A(s).  From  the  dispersion  relation  we  can  obtain  criteria 
for  a post-flare  loop  to  be  susceptible  to  thermal  instability. 


2 


const 


(1.9^ 


to  first  order: 


The  dispersion  relation  is  obtained  from  the  requirement  that  the 
determinant  of  the  coefficient  matrix  in  equation  (1.10)  vanishes. 

As  pointed  out  by  Field,  the  dispersion  relation  has  three  solutions: 
one  mode  corresponding  to  isobaric  condensations  and  the  other  tvo 
corresponding  to  the  growth  of  sound  waves.  In  this  paper  we  will 
only  be  interested  in  the  growth  of  condensations.  We  note  from 
equation  (1.3)  Chat  if  the  speed  of  sound  is  large  so  that  the  acceleration 
term  is  negligible,  then  we  can  replace  the  term  <)1Pq)  in  the  second 
row  of  equation  (1.10)  by  zero.  This  corresponds  to  the  case  of 
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Isobaric  condensations.  Therefore,  to  Isolate  cases,  we  expand  the 


determinant  of  the  coefficient  matrix  about  the  second  row,  yielding: 


o 

- ri  P 


Uy-i)t  p / (v-i)T  '‘1  ^ 


T 


+ Kk 


-)=  0 , (1. 


13) 


where  we  have  dropped  the  subscript  "0".  The  first  term  in  (1.13)  is 
the  contribution  due  to  the  acceleration  term  in  the  force  equation. 
If  this  vanishes  then  the  dispersion  relation  reduces  to 


(1.14) 


For  instability,  we  need  the  real  part  of  T)  to  be  positive.  From 
equations  (1.11)  and  (l.l4),  RgCfl)  >0  implies 

Xj  - I . (1.15) 


Equation  (1.15)  is  identical  to  the  one  derived  by  Field  (I965).  The 
addition  of  a variable  area  has  no  effect  on  the  growth  of  condensations 
other  than  to  add  an  oscillatory  factor  to  the  time  dependence.  If 
conduction  is  negligible,  then  equation  (1.15)  reduces  to 


dT  ■ T dP 


0 


(1.16) 


This  relation  was  originally  derived  by  Weymann  (I960).  Note  that 
Parker's  simpler  condition,  Xp  < 0,  implies  that  the  radiative  loss 
function  for  free-free  emission,  equation  (1.1),  is  stable;  whereas, 
under  condition  (I.I6)  it  is  unstable.  However,  we  find  that  in  the 
corona,  whenever  the  temperature  is  high  enough  so  that  bremsstrahlung 
dominates  the  radiation  losses,  the  conductivity  is  usually  so  large 
that  condition  (I.I5)  implies  stability. 
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We  note  from  equation  (1.1^)  that  the  effect  of  conduction  is  to 
increase  the  thermal  stability  of  the  plasma.  This  is  to  be  expected, 
since  conduction  always  acts  to  reduce  temperature  differences.  For 
given  values  of  n and  T,  equation  (1.1^)  Implies  that  there  exists  a 
maximum  wavenumber,  k^,  such  that  perturbations  of  larger  wavenumber 
are  damped  by  conduction.  Letting 


X = 2n/k 
c c 


and  using  the  Cox  and  Tucker  (I969)  form  for  the  cooling  function. 


(1.17) 


J = n'^A(T)  , 


(1.18) 


where  A “ T , and  a is  of  order  unity,  we  obtain  for  X^, 

' (1.19) 

V 3n  A 

Hence,  one  of  the  criteria  that  a perturbation  must  satisfy  in  order 
for  it  to  grow  is  that 


X ^ X . 
c 


Equation  (I.I9)  provides  a lower  limit  on  the  wavelength  of  an 


(1.20) 


unstable  perturbation.  Field  noted  that  the  dispersion  relation  (1.13) 
also  implies  an  upper  limit  on  the  wavelength  of  an  unstable  density 
perturbation.  This  is  because  the  growth  rate  of  a density  perturbation 
cannot  be  faster  than  the  speed  of  sound  travel  time  across  the  pertur- 
bation; 


T - 1/kc  . 
s 


(1.21) 
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Therefore,  flnlte-speed-o£  sound  effects  will  decrease  the  growth  rate 
if  the  cooling  time, 

(1.22) 

is  less  than  T^.  This  implies  another  criterion  on  the  wavelength  of 
an  unstable  perturbation, 

X . (1.23) 

n“^A 

Note,  however,  that  this  restriction  applies  only  to  the  growth  of 
density  perturbations.  Temperature  perturbations  of  arbitrary  wave- 
length will  grow  on  time  scales  of  the  cooling  time  as  long  as  condition 
(1.20)  is  satisfied.  In  fact,  we  calculate  in  Section  that  temperature 
variations  in  the  post-flare  corona  grow  at  a much  faster  rate  than  density 
variations  due  to  the  finite  speed  of  sound. 
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2.  POST-FLARE  REGIONS 
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2*1  Solar  Prominence ■ 

Along  with  aunapots  and  flarea,  prominencea  are  one  of  the  most 
cooBton  manlfaa tat Iona  of  aolar  activity.  In  general,  a prominence 
la  almply  a region  of  cold,  dense  plasma  that  Is  visible  In  Ha  and  Is 
located  at  least  partially  In  the  corona.  Since  the  corona  usually 
has  a temperature  In  excess  of  10  K,  there  Is  over  two  orders  of 
magnitude  temperature  difference  between  a prominence  and  Its 
surroundings. 

Prominences  have  a very  wide  range  of  morphology;  however,  it  Is 
possible  to  separate  them  Into  two  categories:  quiescent  prominences 
which  are  large  structures  (size  scales  up  to  10^^  cm),  that 
appear  stable  for  months  and  occur  In  quiet  regions;  and  active  promi- 
nences which  are  short-lived  events  (lifetimes  less  than  hours),  and 
are  believed  to  be  flare  related.  (A  thorough  and  detailed  description 
of  the  observed  properties  of  solar  prominences  Is  given  by  Tandberg- 
Hanssen  (197^)  In  his  book  on  the  subject.)  In  this  paper  we  will  be 
Interested  mainly  In  active  loop  prominences.  We  Include  coronal  rain 
and  knots  (Tandberg-Hanssen,  1974,  p.  3)  ••  part  of  this  phenomenon. 

Dodson  (1961)  first  proposed  that  loop  prominence  systems  are 
flare  related  events.  A careful  Investigation  of  the  association 
between  loop  prominence  systems  and  flarea  waa  made  by  Brueek  (1964a). 

Ha  concludes  that  they  occur  only  after  great  flares  and  especially  after 
flares  that  emit  energetic  protons  (Bruzek,  1964b;  Svestka,  I968). 

Teske  (1971)  found  that  loop  prominence  systems  radiate  In  soft  x-rays 
Indicating  that,  In  common  with  post-flare  regions,  they  contain  a 
high  temperature  plasma,  T ^ 10^  K.  A similar  observation  was  noted 
by  waldmeler  (1973).  9 


We  believe  that  loop  prominences  are  a direct  raault  of  thermal 
Instability  In  hot,  danse  post-flare  regions.  They  occur  preferentially 
In  these  regions  of  the  corona  because  the  coronal  density  Is  greatly 
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enhanced  by  a flare.  The  Important  features  of  flare  regions  and  of 
loop  prominences  are  briefly  described  In  the  following  sections.  A 
detailed  description  of  the  flare  corona.  Including  some  of  the  recent 
Skylab  data  Is  available  In  the  book  by  Svestka  (1976). 

2.2  Description  of  the  Post-Flare  Corona 

2.2.1  Magnetic  Field 

The  geometric  structure  of  prominences  and  other  features  In  the 
corona  is  almost  wholly  determined  by  the  magnetic  field.  This  Is 
because  the  field  is  strong  enough  to  prevent  any  mass  or  heat  flow 
perpendicular  to  field  lines  (Section  Hence,  temperature  and 

density  gradients  perpendicular  to  the  field  tend  to  be  much  larger 
than  those  parallel.  In  addition,  mass  motion  will  be  directed  along 
field  lines. 

The  observed  geometry  of  the  field  Is  a loop  or  a series  of  loops; 
each  loop  Is  believed  to  represent  a magnetic  flux  tube.  This  Is  true 
If  a flare  region  Is  observed  either  In  soft  x-rays  (Kahler  et  al., 

197^;  Vorpahl  et  al.,  1975),  or  in  HOT  as  a loop  prominence  system. 

The  magnetic  field  strength  In  loop  prominences  appears  to  be  large, 
as  would  be  expected  for  flare  regions.  Zlrln  (I96I)  measured  fields 
of  up  to  200  G,  and  Hyder  (196^)  obtained  a value  of  60  - 30  G for 
the  field  In  a loop  prominence.  Similar  values  were  observed  by  Roy 
(1972)  and  Rust  and  Bar  (1973)*  Tandberg-Hanssen  (197^,  p.  43)  estimates 
that  the  typical  field  In  a loop  prominence  Is  from  20  - 100  G. 
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Several  observers  have  compared  the  coronal  field  lines  with  those 
of  a potential  field  computed  from  observations  of  the  photospheric 
field.  Rust  and  Roy  (1971).  Roy  (1972)  and  Rust  and  Bar  (1973).  They 
find  that  the  field  in  the  post-flare  corona  can  be  well  approximated 
as  being  potential,  especially  at  large  heights.  This  result  is  in 
agreement  with  the  observation  that  the  field  in  loop  prominence  systems 
is  static.  Bruzek  (1964a)  finds  that  individual  loops  do  not  alter 
their  shape  during  the  prominence  event. 

In  summary,  we  find  that  the  magnetic  field  in  the  post-flare 
corona  is  large,  ~ ^0  G,  and  appears  to  be  current-free. 


2.2.2  Temperature 

Observations  (Thomas  and  Teske,  1971)  imply  that  all  flares  emit 
soft  x-rays.  The  spectra  of  the  soft  x-ray  bursts  Indicate  that  they 

7 

are  due  to  bremsst rah lung  from  a thermal  plasma  at  temperatures  ^10  K. 
Typically  these  bursts  have  rise  times  of  order  minutes  and  decay  times 
of  order  tens  of  minutes.  The  peak  in  temperature  is  usually  observed 
to  occur  before  the  peak  in  flux  (Horan,  1971;  Kahler  et  al.,  1970), 
which  implies  that  the  density  of  the  soft  x-ray  emitting  plasma  is 
rising  while  the  temperature  drops. 

The  dominant  cooling  mechanism  for  the  hot  coronal  plasma  is  believed 
to  be  conduction  to  the  chromosphere,  Culhane  et  al.  (1970),  Moore  and 
Datlowe  (1973)*  A related  process  is  the  evaporation  of  chromospheric 


material.  Although  this  has  not  been  seen  directly,  it  has  been  proposed 
by  several  observers  as  the  mechanism  by  which  the  density  in  the  flare 
corona  is  enhanced  (Neupert,  I968;  Hudson  and  Ohkl,  1972;  and  Neupert 
et  al.,  1974). 
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The  value  of  the  temperature  in  flare  regions  is  generally  estimated 
either  from  the  continuum  or  from  lines  In  the  soft  x-ray  spectrum. 

These  measurements  Indicate  that  the  maximum  temperature  Is  usually 
(2-4)  X 10'  K,  but  that  the  flare  corona  Is  highly  Inhomogeneous  (Neupert 
et  al.,  1974;  Dere  et  al.,  1974).  This  Is  not  surprising  since,  if  a 
flare  region  consists  of  several  loops,  adjacent  loops  are  thermally 
insulated  from  one  another  by  the  magnetic  field.  In  fact,  observations 
indicate  that  loops  of  temperature  > 10^  K can  be  seen  in  the  immediate 
vicinity  of  Ha  loops,  T~  10^  K,  (Fisher,  I97I,  1974;  McCabe,  1973). 

If  conduction  is  indeed  the  dominant  cooling  process,  then  we  can  see 
that  large  temperature  differences  would  naturally  arise.  Even  if  all 
the  loops  in  a flare  region  were  heated  to  the  same  initial  temperature, 
small  loops  would  cool  before  large  ones  since  the  conductive  cooling 
rate  varies  as  L , where  L is  the  loop  length  (Section  ^.2.2). 

2.2.3  Density 

The  density  in  soft  x-ray  emitting  regions  is  not  as  accurately 
known  as  the  temperature.  From  the  observation  of  certain  lines  of 

highly  ionized  elements,  e.g.  Fe  XXV,  one  can  safely  conclude  that 

’ 7 

temperatures  of  at  least  10  K are  present;  however,  there  does  not 
seem  to  be  any  model  independent  method  of  obtaining  values  for  the 
density.  Part  of  the  difficulty  is  that,  like  the  temperature,  there 
is  no  reason  to  expect  a unique  density  for  a post-flare  region. 

From  estimates  of  the  emission  measure,  values  for  the  density  of 
10^^  - 10^^  cm"^  are  usually  given,  e.g.  Cheng  and  Widing  (1975)« 

However,  Craig  and  Brown  (1976)  point  out  that  temperature  inhomo- 
geneities may  lead  to  large  errors  in  the  deduced  magnitudes  of  both 
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density  and  temperature.  In  general,  one  may  say  that  from  soft  x-ray 

observations,  n Is  at  least  10^*^  cm  ^ and  may  be  as  high  as  10^^  cm  ^ 

In  flare  coronal  regions  (Svestka,  1976,  p.  138). 

In  loop  prominences  visible  In  H(V  the  density  may  be  more  easily 

obtained  because  the  plasma  Is  usually  optically  thick  In  the  hydrogen 

lines.  From  the  Stark  broadening  of  the  Balmer  lines,  Hlrayama  (1972) 
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obtains  densities  ranging  from  10  - 10  cm  . Jefferies  and  Orrall 

(1961)  use  the  scattering  of  photospherlc  light  to  calculate  the  electron 

11  -3 

density.  They  obtain  a value  of  2 X 10  cm  In  a loop  prominence. 

In  summary,  we  conclude  that  the  density  In  loops  Is  measured  to 

11  -3 

be  ^10  cm  , while  soft  x-ray  observations  Indicate  that  the  density 
In  a general  post-flare  region  Is  ^ 10^*^  cm  Since  loop  prominence 
systems  are  only  associated  with  large  flares.  It  may  well  be  that  the 
higher  value  Is  more  accurate  for  these  particular  flare  regions. 

2.2.U  Morphology  of  Loop  Prominences 

The  evolution  of  a loop  prominence  system  has  been  studied  In  detail 

by  Bruzek  (1964a).  He  finds  that  these  systems  first  appear  In  Hnr  several 

g 

minutes  after  flare  onset  as  a low  mound,  heights  < 10  cm.  The  system 
then  evolves  as  a series  of  loops  appearing  at  successively  higher 
heights  reaching  up  to  I.5  X 10^*^  cm.  The  apparent  rate  of  expansion  of 
the  system  Is  of  order  a few  km/ sec;  however.  Individual  loops  do  not 
expand.  The  lifetime  of  a typical  loop  Is  approximately  30  minutes  and 
the  complete  system  may  last  well  over  12  hours.  Bruzek  has  observed 
certain  systems  lasting  up  to  3 successive  days. 

Individual  loops  generally  appear  to  begin  as  a bright  small  knot 
at  the  top  of  the  loop,  although  knots  can  sometimes  form  at  the  sides 
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as  well  (Bruzek  and  Kuperus,  1972).  The  knot  then  expands  into  a loop 
with  material  streaming  down  both  legs.  This  part  of  the  evolution  is 
very  similar  to  "coronal  rain"  (Tandberg-Hanssen,  197^»  p.  29).  The 
velocity  of  the  material  at  the  loop  base  is  of  order  10  km/sec,  in 
agreement  with  free-fall  velocities.  The  lifetime  of  a knot  is  of 
order  10  minutes;  Bruzek  and  Kuperus  (1972)  quote  3 ■ 10  minutes,  while 
Tandberg-Hanssen  (197^ » p.  5)  states  I5  minutes. 

2.3  Theories  of  Loop  Prominence  Systems 
2.3.1  Compression-Condensation  Model 

Although  Parker  (1953)  originally  proposed  thermal  instability  as 
a mechanism  for  the  formation  of  quiescent  prominences,  he  discounted  it 
as  a possible  explanation  for  active  prominences  because  their  lifetime 
is  much  smaller  than  the  growth  time  of  thermal  perturbations  at  coronal 

Q 

densities,  < 10  cm  . In  addition,  the  size  scale  implied  by  equation 
(1.20)  for  coronal  densities  is  much  larger  than  active  prominence  size 
scales.  Field's  (I965)  more  exact  analysis  confirmed  this  result. 

Kleczek  (195^)  pointed  out  that  if  the  coronal  plasma  were  sufficiently 
compressed  by  magnetic  fields,  then  the  radiative  cooling  time  may  be 
decreased  enough  to  agree  with  loop  prominence  time  scales.  Hence,  a 

9 "3 

model  based  on  an  initial  strong  compression  (from  n » 10^  cm  to 
n es  10  cm  has  been  investigated  by  several  authors:  Kleczek  (I958), 
Lust  and  Zirln  (i960),  Shklovskli  (I965),  Olson  and  Lykoudis  (I967), 
and  De  ^1973). 

Several  objections  can  be  raised  to  this  model.  Jefferies  and 
Orrall  (1965b)  pointed  out  that  the  amount  of  matter  in  a loop  prominence 

16 

system,  ~ 10*^  gms,  is  too  large  to  condense  out  of  the  corona,  even  in 
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a coronal  condensation.  Olson  and  Lykoudis  (I967),  themselves  agree 

9 

that  the  compression  model  cannot  explain  loops  larger  than  ~ 3 X 10  cm, 
while  loops  larger  than  10^^  cm  have  often  been  observed  (Bruzek,  1964a). 

Teske  (1971)  argues  that  the  compression  model  is  not  consistent  with 
his  observations  of  soft  x-rays  from  loop  prominence  systems.  Finally, 
this  model  requires  large  changes  in  the  magnetic  field  in  order  to 
compress  the  plama.  However,  observations  indicate  that  the  field  is 
static  during  a loop  prominence  event,  in  fact,  several  authors  find 
that  it  is  potential  (Section  2.2.1).  Therefore,  we  conclude  from 
these  arguments  that  the  compression-condensation  model  is  not  a likely 
explanation  for  the  formation  of  loop  prominences. 

2.3.2  Fast  Particle  Model 

" 1 

Jefferies  and  Orrall  (1965a,  1965b)  have  proposed  a model  for  loop  i 

1 

1 

prominence  systems  in  which  the  prominence  material  is  supplied  by  fast 

Q J 

protons,  Vp  ~ 10°  cm/ sec.  They  argue  that  the  large  mass  required  for 
a loop  system  is  stored  as  10  keV  protons  on  large  loops  extending  into  . 

the  outer  corona.  These  particles  are  supposed  to  migrate  down  to  , 

i 

lower  field  lines  where  they  thermallze  near  the  top  of  a loop  to  produce  | 

the  bright  knots  that  are  observed.  | 

The  major  objection  to  this  proposal  is  that  it  is  difficult  to  | 

imagine  how  the  protons  can  migrate  across  the  field,  and  in  such  a | 

way  as  to  reproduce  the  evolution  of  a loop  prominence  system.  Particles  | 

stored  in  the  outer  corona  would  have  to  migrate  first  down  to  the 
lowest  loops  and  then  to  higher  loops.  Unless  some  mechanism  whereby 
this  may  be  done  is  proposed,  the  fast  particle  model  cannot  be  truly 
evaluated. 
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Goldsmith  (1971)  and  Sturrock  (1973)  have  proposed  that  the  mass 

for  a loop  prominence  system  is  due  to  chromospheric  evaporation  by 

flare  heating.  They  argue  that  since  soft  x-ray  observations  indicate 

that  the  density  in  post-flare  regions  is  greatly  enhanced  over  usual 

coronal  values,  a loop  prominence  system  may  well  be  simply  the  final 

stage  of  the  evolution  of  the  dense  flare  corona.  The  appearance  of 

loops  at  successively  greater  heights  may  be  explained  as  due  to  the 

decrease  in  the  conductive  cooling  rate  with  height. 

The  idea  of  chromospheric  evaporation  has  another  attractive  feature 

in  that  one  does  not  have  to  invoke  magnetic  field  compression  or  other 

exotic  mechanisms  in  order  to  obtain  large  enough  densities  for  thermal 

instability  to  be  effective.  Goldsmith  has  calculated  the  response  of 

a coronal  loop  to  an  Isobarlc  perturbation  in  its  Initial  temperature 

6 1 1 • 

and  density  profile.  He  finds  that  for  Tq  = 5 X 10  K,  and  n^  = 10  cm 

, 9 

a 5%  amplitude  perturbation  of  wavelength,  3 X 10  cm  does  grow  to 

produce  temperature  differences  of  over  an  order  of  magnitude.  However, 
these  large  temperature  differences  last  for  only  < 100  sec,  which  is 
small  compared  to  the  lifetime  of  knots,  ~ 10  minutes.  In  addition, 
his  calculations  do  not  include  effects  due  to  gravity,  conduction  to 
the  chromosphere,  or  variation  of  cross-sectional  area  of  the  loop. 

7 

Therefore,  the  Initial  cooling  of  the  10  K plasma,  and  its  effect  on 

the  formation  of  condensations  cannot  be  investigated  and  the  dynamics 
k 

of  the  cool,  10  K,  plasma  cannot  be  calculated.  Hence,  this  model  cannot 
be  directly  compared  with  the  observations  of  loop  prominences. 

We  believe  that  the  general  features  of  the  evaporation  model 
proposed  by  these  authors  are  correct.  In  this  model  loop  prominences 
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are  a natural  part  of  the  cooling  process  and  are  not  due  to  some 
unobserved  phenomena.  Therefore,  in  order  to  critically  test  the 
model,  we  calculate  the  evolution  of  a post-flare  loop  including  all 
the  effects  described  above;  and  compare  our  results  with  observations. 
The  model  of  a loop  prominence  that  we  use  for  our  calculation  is 
described  in  detail  in  the  following  chapter. 
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3.  MODEL  OF  A LOOP  PROMINENCE 

3.1  Introduction 

In  this  chapter  we  describe  a one-dimensional  (only  variations 
parallel  to  B are  considered),  single-fluid  model  of  a loop  prominence 

1 

which  we  have  analyzed  numerically.  The  physical  assumptions  that  have 
been  made  and  their  justification  are  discussed  below. 

3.2  Single-Fluid  Assumption 

I 

The  temperatures  and  densities  used  in  the  model  are  those  observed  j 

for  loop  systems,  j 

j 

3 X 10^  K T ^ lo"^  K , (3.1)  I 

11  -3 

and  n ~ 10  cm  . We  assume  that  the  plasma  can  be  described  by  single- 

fluid  equations.  This  approximation  is  valid  only  if  the  characteristic  ‘ 

size  scales  and  time  scales  of  the  fluid  are  larger  than  the  mean  free 

path  and  the  relaxation  time  respectively.  Spitzer  (I962)  has  calculated  j 

the  mean  free  path  of  electrons  and  protons  in  a fully  ionized  hydrogen  1 

plasma  and  the  electron-proton  relaxation  time, 

Xj  ««  10^  T^  n"^  cm.  , (3.2)  : 

and 

T (^20  T^^^  n"^  sec.  (3.3) 

The  length  and  time  scales  over  which  significant  changes  can  occur  in 

a loop  are  the  thermal  instability  size  scale,  X^,  equation  (I.I9)  and 

the  cooling  time  equation  (1.22).  For  the  values  of  T in  (3*1)»  i 

we  find  that  X,/X  and  T /t  < < 1.  ] 

t c ep  c i 
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3.3  Radiative  Losses 

Since  hydrogen  comprises  ~ 90^  of  the  ions  in  the  coronal  plasma, 
the  effects  of  the  heavier  trace  elements  are  negligible  except  for 
their  contribution  to  the  radiation  loss  rate  at  temperatures  above 
~ 3 X 10^  K.  At  high  temperatures,  10^  K),  the  radiation  cooling 
is  due  almost  entirely  to  lines  in  heavy  ions  such  as  0^.  These  losses 
can  be  represented  as  a heat  sink  term  in  the  energy  equation,  thus,  we 
approximate  the  loop  material  by  the  fluid  equations  for  a pure  hydrogen 
plasma  but  with  a cooling  function,  X,  which  includes  the  radiative 
losses  of  the  heavy  elements. 

The  form  of  X used  is  that  derived  by  Cox  and  Tucker  (1969).  In 
their  calculations  they  assumed  that  the  plasma  is  of  cosmic  abundances, 
is  optically  thin  to  all  radiation,  and  is  in  statistical  equilibrium 
with  the  states  of  ionization  determined  by  a balance  of  radiative  and 
dielectric  recombination  with  collisional  ionization.  All  the  radiative 
and  ionization  processes  considered  vary  as  the  electron  density  and 
ion  density,  hence  the  cooling  function  can  be  written  as: 

X(n,T)  - n^  nA(T)  erg  cm  ^ s ^ , (3.4) 

where  A(T)  has  been  tabulated  by  Cox  and  Tucker,  Figure  1.1. 

4 

The  assumptions  above  are  valid  for  loop  plasma  above  ~ 3 X 10  K 
where  hydrogen  is  almost  fully  ionized,  but  they  break  down  at  temper- 
atures below  this.  Once  a significant  amount  of  neutral  hydrogen 
forms,  the  plasma  will  not  be  optically  thin  to  hydrogen  lines.  The 
mean  free  path  of  a line  photon  is  given  by, 

X = 1/tij^  a , (3.5) 
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where  a Is  the  cross-section  for  absorption.  For  Ly  Of,  a ~ 10  cm  , 

3 11  -3 

hence  the  mean  free  path  is  ~ 10  cm  at  densities  of  10  cm  . This 
means  that  if  a slab  of  loop  plasma  only  1 km.  thick  were  to  cool  to 

4 2 

temperatures  ~ 10  K,  the  optical  depth  of  the  slab  would  be  ^ 10  . 
Clearly,  condensation  in  loops  cannot  be  optically  thin  in  at  least 
Ly  Of,  hence  , in  order  to  follow  the  cooling  at  low  temperatures  in 
detail  one  must  Include  radiation  transfer  effects. 

Another  assumption  that  breaks  down  at  low  temperatures  is  that 
the  state  of  ionization  is  in  equilibrium.  Consider  a plasma  that  is 
cooling  from  a very  high  temperature,  > 10^  K,  so  that  initially  all 
the  hydrogen  is  ionized.  Assuming  that  the  Cox  and  Tucker  conditions 
are  applicable,  it  will  eventually  cool  to  1.6  X 10^  K with  percentage 
of  ionization,  n^/n  ss  ^0%.  Now  in  order  to  Justify  using  equilibrium 
equations  to  determine  the  ionization  state,  we  must  have  the  cooling 
time  much  greater  than  the  recombination  time,  otherwise  the  temperature 
will  change  before  equilibrium  can  be  established.  However,  with  the 
cooling  rate  and  the  recombination  coefficient  given  by  Cox  and  Tucker, 


we  find  that  at  T ^ 1.6  X 10^  K, 


Vc  ■ (^)  ' • 


(3.6) 


hence  the  assumption  of  statistical  equilibrium  is  invalid  at  these 


temperatures. 


The  arguments  above  imply  that  in  order  to  calculate  the  behavior 

If, 

of  the  plasma  below  the  3 X 10  K level,  we  would  need  to  follow  it  on 
very  small  size  and  time  scales  and  include  radiation  transfer.  There- 


fore, mainly  for  ease  of  calculation,  we  assume  in  our  model  that  the 


ll 

temperature  does  not  fall  below  3 X 10  K.  This  assumption  is  not 

unjustified  since  we  expect  that  for  the  physical  conditions  present 

in  loops,  the  Important  heat  loss  mechanisms,  conduction  and  radiation, 

4 

essentially  vanish  by  ~ 10  K.  Radiation  losses  become  small  since 
the  temperature  is  ^ 1/10  the  threshold  energy  required  to  excite 
hydrogen  from  its  ground  state,  and  the  number  of  free  electrons  is 
reduced  due  to  recombinations.  Also,  optical  depth  effects  result 
in  a decrease  of  emission.  Conduction  becomes  negligible  at  low 
temperatures  because  the  coefficient  of  conduction  for  a fully  ionized 
plasma  varies  as  a high  power  of  T (Spitzer,  I962): 

K = 10”^  rj.5/2  gjgg  cm”^  k"^  sec"^  . (3.7) 

By  fixing  a minimum  temperature  we  do  ignore  the  profile  of  the 

4 4 

plasma  between  10  - 3 X 10  K.  However,  we  are  primarily  interested 

in  determining  whether  the  thermal  instability  mechanism  can  produce 
condensations  in  loops.  As  far  as  the  hot  loop  plasma,  > 10^  K,  is  con- 

4 4 

earned,  the  difference  between  10  and  3 X 10  K is  negligible.  The 
main  error  is  that  we  effectively  overestimate  the  pressure  of  a 
condensation,  hence  we  can  say  that  if  condensations  do  form  under 
our  assumptions,  they  will  certainly  form  in  a more  exact  calculation 
and  will  be  cooler  and  smaller. 

3.4  MHD  Equations 

4 

Since  we  restrict  the  loop  temperature  to  T ^ 3 X 10  K,  the  state 
of  ionization  n^/n  > 99'^>  ve  can  take  the  plasma  to  be  fully  ionized. 
The  single-fluid  equations  for  a fully  ionized  hydrogen  plasma  in  the 
presence  of  a gravitational  field  and  a magnetic  field  are  well  known. 
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and  Che  equation  of  state 


The  quantity  c In  (3.10)  Is  the  Internal  energy  per  gram  of  the  plasma 
and  In  our  case. 


me  stress  censor  u and  the  heat  flux  vector  q can  be  broken  up 
Into  two  parts,  one  parallel  and  e perpendicular  to  the  magnetic 
field. 


where  Che  elements  of  W are  terms  like 


For  our  model  the 


rtiti  tg  ilhjiMfitllifti  i" 


components  of  q and  II  perpendicular  to  B are  negligible  compared  to 
those  parallel.  This  is  because  the  diffusion  coefficients  for  a 
particular  direction,  ^ and  ^ depend  on  the  mean  free  path  in  that 
direction.  Parallel  to  the  field  the  relevant  quantity  is  the  collision 
mean  free  path,  and  perpendicular  to  the  field  the  gyroradius  is  important. 
Braginskii  (I965)  obtains. 


T . 


(3.14) 


where  is  the  collision  time  and  lo  is  the  gyrofrequency.  Taking 


B ~ 10^  G,  and  n ~ 10^^  cm  yields. 


(<o  T .f=  10^  B^  n"^  = lO'^^  . 

g col^ 


(3.15) 


For  temperatures  T ^ 10"^  K,  the  ratio  of  parallel  to  perpendicular 

4 

diffusion  is  ^ 10  , hence  we  can  ignore  perpendicular  diffusion.  For 
low  temperatures,  ^ 10^  K,  the  mean  free  path  is  so  short  that  both 
parallel  and  perpendicular  diffusion  are  negligible. 

The  coefficient  of  viscosity,  , is  given  by  Spitzer  (I962) 


= 10  gm  cm  ^ sec  ^ , 


(3.16) 


For  our  values  of  density,  temperature,  and  size  scale,  the  Reynolds 
number  for  the  plasma  is. 


“ % 


(3.17) 


This  large  value  for  means  that  the  viscosity  results  in  a negligible 
dissipation  of  momentum  and  energy,  i.e.  the  plasma  is  approximately 
inviscid.  In  the  numerical  calculations  we  introduce  an  artificial 
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viscosity  in  order  to  take  care  of  the  possible  formation  of  shocks, 
however  we  ignore  the  ordinary  plasma  viscosity. 

3.5  Magnetic  Field  Structure 

In  our  model  it  is  assumed  that  all  quantities  vary  along  B only; 
hence,  one-dimensional  equations  can  be  used.  Justification  for  this 
assumption  is  provided  by  the  following  arguments. 

The  magnetic  field  in  a post-flare  region  is  observed  to  have 


sufficient  strength  to  dominate  the  plasma.  This  may  not  be  true 
immediately  after  a flare;  since,  if  the  plasma  thermal  energy  is 
derived  from  magnetic  energy,  then  probably  they  are  initially  of 

comparable  magnitude.  However,  as  the  plasma  cools  from  initial 

7 7 

temperatures  of  ~ 3 X 10  K to  < 10  K the  magnetic  field  begins  to 
dominate.  For  our  temperatures  and  densities  the  field  strengths 
reguired  are,  B ^ 10  G,  a value  well  within  the  magnitudes  observed 
(Section  2.2.1). 


Since  the  plasma  is  fully  ionized,  the  conductivity  is  high 
enough  to  insure  that  the  "frozen  in"  condition  is  valid,  i.e.  the 


magnetic  flux  through  any  area  moving  with  the  plasma  remains  constant. 
In  our  case  the  plasma  pressure  is  not  sufficient  to  transport  the 


field;  hence,  the  plasma  is  constrained  to  move  only  along  field  lines 
like  beads  on  a rigid  wire.  It  was  shown  in  the  previous  section 
that  the  only  significant  components  of  the  stress  tensor  II  and  the 
heat  flux  q are  also  parallel  to  B.  Therefore,  the  equations  of 
motion  (3*9)  reduce  to  a single  equation  for  mass  flux  along  B,  and 
the  heat  equation  (3.10)  involves  heat  flux  only  along  B as  well. 


i 


►4 


r-1 


In  order  to  obtain  a one -dimens tonal  model  we  must  assume  that 
the  relevant  variables  T,  v and  n do  not  vary  over  some  area,  A, 
perpendicular  to  B.  The  observed  morphology  of  loop  prominences 
(Bruzek,  1964a)  suggests  that  they  form  along  field  lines,  the  loop 
consisting  of  a flux  tube  with  cross-sectional  radius  r^  small 
compared  to  the  loop  length.  For  r^  sufficiently  large  then  the 
arguments  above  imply  that  the  tube  is  effectively  insulated  from  the 
rest  of  the  corona,  since  during  the  lifetime  of  a loop  very  little 
heat  or  mass  is  exchanged  between  the  loop  and  its  surroundings.  This 
will  be  true  when  the  cooling  time  t of  a loop  is  much  less  than  the 
time  scale  for  significant  transverse  conduction,  T^.  The  determining 


factor  is  t/t  , where 
x’ 


-2-  - in-^  r2  -1  ^ 2 

K T ■ 10  ® T " ‘'O  ’ 


(3.18) 
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and  T 10  sec.  from  observations.  Hence  for  B ««  30  G,  T 10  K, 

and  n « 10^^  cm  ~ * 1 for  r^  ■ 10^  cm.  If  the  tube  radius  is  less 
than  10°  cm.  then  conduction  across  the  tube  is  dominant;  and,  the 
temperature  will  be  approximately  constant  over  the  cross-sectional 
area  of  the  tube.  However,  this  means  that  a significant  amount  of 
heat  may  be  exchanged  between  the  tube  and  the  surrounding  corona. 
Conversely,  if  the  radius  r^  is  greater  than  10^  cm.,  then  the  tube  is 
effectively  insulated  from  the  outside  but  significant  temperature 
variations  can  occur  across  the  tube. 

g 

Loop  prominences  have  radii,  r^  10°  cm.,  and  so  we  expect  the 
latter  case  to  hold.  Hence,  there  could  be  small  scale  irregularities 
across  loops,  but  structures  of  this  size  scale,  < 10^  cm.,  would  be 
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very  difficult  to  observe.  In  our  model  we  assume  that  there  are  no 

small  scale  Irregularities;  therefore,  the  loop  is  both  insulated  and 

constant  on  a cross-sectional  area.  This  assumption  allows  us  to  use 

one-dimensional  equations.  Of  course,  there  is  a boundary  between  the 

loop  and  the  surroundings  over  which  large  temperature  differences 

occur;  however,  the  discussion  above  indicates  that  the  width  of  this 

boundary  is  small  compared  to  the  width  of  loop  prominences.  Since 

we  do  not  specify  the  magnitude  of  the  area  of  the  tube  in  our  equations,  1 

the  calculations  also  apply  to  the  case  where  the  tube  area  is  small 

enough  to  insure  that  no  small  scale  irregularities  are  present,  but  j 

the  tube  is  also  Insulated  from  the  surrounding  plasma.  This  situation  | 

might  be  possible  if  the  field  is  non-potential. 

Figure  3*1  illustrates  the  geometry  of  our  model  loop.  We  assume 
a dipole  source  Imbedded  a distance  D below  the  top  of  the  chromosphere, 
and  consider  a flux  tube  of  height  above  the  source  R that  lies  in  a 
plane  perpendicular  to  the  solar  surface.  We  could  easily  consider 
the  tube  to  be  inclined  at  some  angle  from  the  vertical.  This  would 
only  reduce  the  effective  gravitational  acceleration  g by  some  constant 
factor,  cos  (p,  and  would  not  appreciably  affect  our  results.  In  general, 
loop  prominences  are  not  observed  to  occur  at  large  angles  from  the 
vertical. 

Letting  s be  distance  along  the  loop  measured  from  the  top,  we  have 
that  all  variables,  including  the  cross-sectional  area  A,  are  functions 
of  s only.  The  effect  of  the  magnetic  field  on  the  plasma  is  accounted 
for  by  Including  A and  by  permitting  only  heat  and  mass  transfer  along 
s.  Each  flux  tube  is  uniquely  characterized  by  two  parameters:  the 
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size  R and  the  compression  factor  F which  is  defined  as  the  ratio  of 
field  strength  at  the  chromosphere  (s  = ± s,  ) to  that  at  the  top  (s  = O). 


r = BfS|^)/B(s  = 0)  = A(s  = 0)/A(s  = Sj^l 


(3.19) 


3.6  One-Dimensional  Equations 


Modifying  the  MHD  equations  (3.8)  - (3.10)  to  fit  our  one-dimen- 
sional model  and  neglecting  viscosity,  we  obtain  the  following  set  of 
equations ; 


n 


dn  . 1 5 , . X 

AM  ^ ° ’ 


P - ng„ (s ) = 0 , 


(3.20) 


(3.21) 


The  derivation  of  equations  (3.20)  - (3.23),  including  the  viscosity 
term  is  given  in  Appendix  A.  In  these  equations  v is  the  velocity 
along  the  field,  and  g is  the  component  of  the  gravitational  field 
along  B.  The  variable  n is  the  number  density  of  ions  (mostly  protons 
in  our  case);  hence  the  mass  density. 


p = (m^  + m ) n = m^^n 


(3.23) 


the  pressure. 


p = (n  + n ) kT  = 2nkT  , 
e p ' 


(3. 24) 


and  the  internal  energy  per  ion. 


The  area  of  the  tube  is  given  by  A(s);  however,  from  the  equations 


above  we  note  that  the  magnitude  of  A is  unimportant,  only  relative 


values  enter.  Therefore,  without  loss  of  generality  we  set  the  area 


at  the  top  to  unity 


are  sufficient  to  determine  our  four 


Equations 


as  functions  of  t and  s 


variables  v,  t,  n and 


are  known 


functions  of  s;  they  are  completely  determined  by  our  choice  of 


magnetic  field.  In  Appendix  A the  dependence  of  A and  g on  s is 


derived 


where  9 is  the  angle  from  the  vertical  (Figure  3.1  )•  We  obtain 


Although  g and  A are  not  expressed  explicitly  in  terms  of  s,  their 


numerical  value  for  a particular  s can  easily  be  obtained  using 


equations  (3.28)  - (3«30),  Since  the  heights  of  the  loops  we  consider 


are  small  compared  to  a solar  radius,  equation  (3>30)  does  not  include 


volume  of  the  loop  between  the  top  and  the  point  s. 


iJ 

^ 1 


A new  variable,  S(x,t),  must  now  be  Introduced.  It  denotes  the 
Eulerian  position  of  that  slab  of  plasma  labeled  "x".  Converting 
equations  (3.20)  - (3.22)  to  Lagrangean  formulation  and  rearranging, 
we  obtain: 


L-  v(x,t) 

at 


1 ^ P(x,t) 

nijj  A(S)  Sx 


+ g„(s(x,t)) 


(h.3) 


a s (x,t ) 

at 


v(x,t) 


(4.4) 


^ T(x,t)  _ pfx.t ) ^ (a(S)  v(x,t))  n(x.t)  A(x) 

at  " ■ 3k  ax  ' ' ■ 3k 


^ ■ lo'^  a / 2 ^ 

21  • k ax  \ dx  j * 


(4.5) 


with 


p = 2knT 


(4.6) 


and 


(4.7) 


Equations  (4.3)  - (4.7)  are  in  the  form  that  we  use  for  numerical 
calculations.  We  assume  a grid  of  equal  intervals  in  x,  which  is 
equivalent  to  dividing  up  the  loop  into  slabs  of  equal  mass.  The 
behavior  of  the  loop  plasma  is  then  simulated  by  following  the  evolution 
of  these  slabs  as  determined  by  equations  (4.3)  - (4.5). 
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4.2  Boundary  and  Initial  CondltlonB 


t- 


Sr 


In  order  to  determine  the  correct  boundary  conditions  for  the 
system  (4.3)  " (4.7)#  we  first  note  that  there  are  actually  only  two 
distinct  dependent  variables:  the  position  of  the  plasma  S and  the 

temperature  or  internal  energy  T.  Equations  (4.3),  (4.6)  and  (4.7) 
can  be  used  to  eliminate  v,  p and  n,  resulting  in  one  partial  differential 
equation  of  second  order  in  time  for  S and  one  of  first  order  in  time  for 
T.  The  spatial  terms,  however,  are  highly  nonlinear.  The  highest  x 
derivative  for  both  S and  T are  of  second  order.  We  essentially  have 
a hyperbolic  equation  for  S and  a parabolic  equation  for  T in  the  two 


Independent  variables  x and  t.  Suitable  boundary  conditions  for  such 
a system  are  Cauchy  conditions  for  S and  Dirichlet  for  T on  an  open 
region  in  the  x,t  plane.  In  particular,  as  Initial  conditions  we 
specify  T,  S and  ^ on  t = 0,  and  as  spatial  boundary  conditions  we 
specify  T and  S at  both  ends  of  the  loop  for  all  t. 

From  the  geometry  of  the  loop,  Figure  3«1  and  from  observations, 

it  appears  reasonable  to  assume  that  the  loop  is  symmetric  about  the 

vertical  axis.  This  implies  that  T must  be  a symmetric  function  of 

X and  S antisymmetric;  hence  we  need  solve  the  equations  only  for 

X ^x^O.  Atx=0,  i.e.  the  top  of  the  loop,  the  boundary  conditions 
n 

must  be  ^ , q ^^d  S(0,t)  = 0. 

The  conditions  at  the  base  of  the  loop  are  determined  from  physical 
consideration.  We  are  interested  in  the  cooling  phase  of  flares, 
specifically  when  radiative  losses  dominate  the  cooling.  During  this 
phase  we  do  not  expect  to  have  significant  heating  or  evaporation  of 
chromospheric  material;  therefore,  we  take  the  temperature  at  the  base 
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of  our  model  to  be  that  in  the  upper  chromosphere,  T (Sj^)  ««  3 X 10^  K, 
and  to  remain  constant  during  our  calculation.  We  do  Include  the 
conduction  of  heat  from  the  loop  to  the  base;  however,  we  assume  that 
this  energy  can  be  radiated  away  by  the  base  without  appreciably 
affecting  Its  temperature. 

This  situation  may  appear  Inconsistent  since  on  the  one  hand,  we 
have  that  the  cooling  rate  of  the  base  is  high  enough  to  radiate  all 
of  the  energy  conducted  down  from  the  loop  yet,  on  the  other  hand,  the 
rate  Is  so  slow  that  the  base  temperature  does  not  change  during  the 
cooling  time  of  the  loop.  The  resolution  to  this  Inconsistency  Is 
In  the  behavior  of  the  radiative  losses  at  low  temperatures.  As 
discussed  in  Section  3.3»  the  loss  rate  decreases  very  rapidly  for 
low  temperatures.  Assuming  the  Cox  and  Tucker  loss  rate,  we  obtain 
for  n ~ 10^^  cm  ^ and  T ~ 3 X 10^  K 

h _5  _1 

n^nA(T)  ~ 10  ergs  cm  s 

This  radiative  rate  Is  high  enough  to  dissipate  an  Incoming  energy 
9 -2-1 

flux  of  10  ergs  cm  s In  a thickness  of  only  1 km.  But  at 

T ~ 10^  K,  n^nA(T)  < 10  ^ ergs  cm  ^ s which  Implies  a cooling  time, 

Tc  ~ ^ 10^  s.  Therefore,  the  base  plasma  is  "trapped"  between 

4 ® 4 

10  - 3 X 10  K;  the  cooling  rate  at  the  lower  temperature  is  too  low 

to  allow  a significant  decrease,  while  the  losses  at  the  upper  value 
are  too  high  to  permit  an  Increase  In  temperature.  The  exact  profile 

of  the  base  temperature  Is  unimportant  to  our  model;  hence,  for 

4 

convenience  we  take  it  to  be  constant  at  3 X 10  K* 


As  a boundary  condition  on  S,  the  position  of  the  base  is  assumed 

to  be  fixed  during  the  cooling  process.  This  is  a good  approximation 

since  the  gravitational  scale  height  of  the  chromosphere < is  very  small 

4 kT 

compared  to  the  size  of  typical  loops.  At  T ~ 10  K,  ~ 400  km; 

therefore,  the  chromospheric  material  will  move  very  little  in  response 
to  pressure  changes  in  the  loop  plasma. 

The  initial  conditions  for  our  model  loop  are  somewhat  arbitrary. 
Order  of  magnitude  values  of  the  temperature  and  density  are  set  by 

7 11  -3 

observations,  T / K and  n ~ 10  cm  ; however,  the  detailed  profiles 
are  not  known.  As.  !;onditions  on  the  velocity  and  pressure,  we  assume 
that  the  plasma  is  initially  at  rest. 


'(s,  0)  = 0 , 


(4.8) 


and  in  hydrostatic  equilibrium. 


t - 0 


"Hi 


(4.9) 


Neither  of  these  assumptions  is  probably  valid  immediately  after  the 
loop  fills  with  evaporated  chromospheric  material;  however,  we  expect 
these  conditions  to  be  achieved  in  a time  short  compared  to  the  cooling 

3 Q 

time,  2 10  secs.  For  loop  heights  of  ~ 3 X 10^  cm.  the  travel  time 

2 

across  the  loop  at  the  sound  speed  is  only  ~ 10  s;  hence,  plasma  motion 
will  cancel  unbalanced  pressure  gradients  before  appreciable  cooling 
can  take  place.  This  motion  will  in  turn  be  damped  by  thermal  conduction; 

g 

the  dissipation  time  for  wavelengths  < 3 X 10°  cm.  is  < 30  seconds. 

9 

Oscillations  of  long  wavelengths,  10  cm.,  corresponding  to  a funda- 


mental mode  of  the  tube  could  persist;  however,  our  calculations 


(Section  5.U.3)  indicate  that  they  do  not  significantly  affect  the 
formation  of  condensations. 
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The  important  parameters  of  the  initial  temperature  profile  are 
the  average  temperature,  which  is  determined  from  soft  x-ray  observations, 
and  the  heat  flux  to  the  chromosphere.  Only  an  upper  bound  can  be 
placed  on  the  latter;  the  downward  heat  flux  should  not  be  more  than  the 
total  radiative  rate  during  maximum  flare  brightness.  (This  argument, 
of  course,  assumes  that  a significant  portion  of  the  flux  is  not  dissi- 
pated by  evaporation.)  Since  flares  are  not  usually  seen  in  white  light, 
we  can  place  an  upper  bound  of  ~ 10^^  ergs  cm  ^ s ^ on  the  heat  flux 
into  the  chromosphere;  however,  our  calculations  (Section  ^.2),  indicate 
that  the  flux  is  usually  much  less  than  this  value. 

Assuming  that  the  heat  flux  is  of  the  ^ame  order  as  the  chromospheric 
7 -2-1 

Hot  flux,  F ~ 10  ergs  cm  s , we  obtain  for  the  temperature  scale  height; 


(^4.10) 
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hence  = 3 10  cm.  at  T » 10  K.  Since  this  scale  height  is  much 

laiger  than  the  size  of  a loop,  only  a very  slight  temperature  difference 
over  the  loop  length  is  needed  to  carry  the  flux.  Also,  the  scale  height 
has  a strong  dependence  on  T.  For  low  temperatures,  T ^ lo'^  K,  equation 

g 

(4.10)  yields  < 10  cm.,  which  is  much  less  than  typical  loop  lengths. 
From  the  arguments  above,  we  expect  that  most  of  the  loop  volume  will 
be  at  ~ 10  K,  with  only  a small  fraction  of  the  plasma  at  lower  temper- 
ature. Therefore,  we  approximate  the  profile  by  assuming  that  initially 

7 

the  temperature  is  uniform  at  ~ 10  K,  i.e. 

T ( 8,  0)  * const.  . (4.11) 
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Equation  (4.9)  can  now  be  used  to  determine  the  pressure  profile. 


p(s,  0)  - p 


b f-^  / 

L n 


g„  (•)  <*»  . 


(4.12) 


Howeyer,  for  T ~ 10'  K,  the  gravitational  scale  height 


H - 10^^  cm.  , 

* “h*ii 


(4.13) 


which  is  much  larger  than  the  loop  length.  Hence,  equation  (4.12)  can 
be  well  approximated  by 


p(s,  0)  - p - const. 


(4.14) 


Equation  (4.11)  and  (4.l4)  imply  that 


n(s,  O)  « const. 


(4.15) 


Using  the  Initial  density  profile  (4.15)  and  equation  (4.1),  we 
can  express  the  initial  conditions  (4.8),  (4.11)  and  (4.15)  in  terms  of 
X.  In  this  case,  the  procedure  is  trivial;  all  three  initial  profiles 
are  uniform  in  x.  They  have  the  obvious  advantage  of  being  very  simple. 
Perturbation  can  easily  be  introduced  into  these  profiles.  There  is 

the  difficulty,  however,  of  matching  them  with  the  boundary  conditions 

4 7 

at  the  base.  Since  the  base  is  at  3 X 10  K while  the  loop  is  at  10 
we,  essentially,  have  replaced  the  steep  gradient  in  T by  a discontinuity. 
Therefore,  we  must  modify  the  boundary  conditions  at  the  base. 

The  Important  quantity  is  the  heat  flux  out  the  loop.  We  use  the 
following  method  to  determine  it.  Assume  that  the  loop  has  been  divided 
into  a number  of  equal  mass  intervals.  The  position  and  the  temperature 
T|^  of  the  base  are  known  and  assumed  to  be  constant  in  time;  also,  the 
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mean  position  and  temperature  of  the  plasma  in  the  lowest  slab 


$£  and  are  known  at  some  time  t.  The  outward  heat  flux  Is  then 
calculated  by  assuming  that  the  flux  between  and  Is  Independent  of 
X,  l.e.  we  assume  that  for  fixed  t the  heat  flux  across  the  final  slab 
is  a constant  heat  flux  determined  by  its  temperature  and  width.  Thus, 

.-6 


10 


/ 7/2  _ 7/2\/  /^b 

^ ^b  .(s)j 


(4.16) 


where  we  have  used  the  Spitzer  conductivity.  If  the  radiative  losses 

of  the  plasma  between  S,  and  S.  are  small,  then  the  expression  (4.l6) 

r D 

is  valid.  If  the  losses  are  large,  then  we  are  overestimating  the 
outward  heat  flux  since  a significant  portion  of  it  will  be  radiated 
away  before  reaching  the  base.  In  our  calculations,  is  usually 

small  enough  so  that  the  first  situation  is  realized.  In  any  case,  the 
maximum  error  is  only  of  the  order  of  a factor  of  2 in  the  total  energy 
losses  of  the  last  slab.  Also,  our  loss  rate  function  A(T)  is  probably 
only  of  that  order  of  accuracy.  Since  we  are  primarily  interested  in 
the  behavior  of  the  plasma  at  the  top  of  a loop,  where  condensations 
form,  we  do  not  expect  that  the  use  of  equation  (4.16)  introduces  a 
significant  error. 

In  summary,  we  have  that  equations  (4.3)  - (4.7)  to  be  solved 
on  the  open  region. 


0 i X < X , and  t ^ 0 , 

m 


with  given  initial  conditions 

T(x,  0)  , S(x,  0)  and  v(x,  0))  , 


(4.18) 
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and  spatial  boundary  conditions 

^T(0,  t)  - 0 , T(x  t)  - 3 X 10^  K , 
dx  “ 

S(0,  t)  » 0 and  S(x  , t)  » S.  ■ const. 

n D 

Note  that  equation  (4.4)  and  the  conditions  above  Imply  that 


(4.19) 


v(0,  t)  - v(x^,  t)  - 0 . 


(4.20) 


4.3  Shocks 

A sticky  problem  In  the  numerical  solution  of  hydrodynamic  equations 
such  as  (4.3)  - (4.^)  is  the  handling  of  shocks;  no  fast  method  for 
treating  them  exactly  exists.  We  use  the  standard  procedure  (e.  g. 
Goldsmith,  1970)  of  introducing  the  van  Neumann-Rlchtmeyer  artificial 
viscosity  term  in  the  equations.  This  term  is  a pseudo-viscous  pressure. 
In  equations  (4.3)  and  (4.5)  the  pressure  p is  replaced  by  (p  + q)  where 
q is  defined  as 


5v 

dx 


< 0 


0 


(4.21) 


where  Ax  is  the  grid  spacing  and  a is  a dimensionless  constant,  a ~ I.5  - 
2,  which  determines  the  thickness  of  the  shock. 

The  effect  of  introducing  this  term  in  the  fluid  equations  has  been 

Investigated  In  detail  (Richtmeyer  and  Morton,  1%7).  A true  shock 

< 

2 

shows  up  as  a Jump  in  the  values  of  n,  T and  v with  a width  of  ~ a 
intervals  and  with  approximately  correct  values  for  the  size  of  the 
jump  and  the  speed  of  the  shock  through  the  fluid.  Because  ^ ^ j Is 
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small  outside  of  the  shock,  q does  not  affect  the  behavior  of  the  plasma 
away  from  the  shock.  In  addition,  since  q Is  defined  to  vanish  for 
positive  velocity  gradients.  It  has  no  effect  on  free  expansion.  Tests 
on  the  artificial  viscosity  by  various  authors  (Rlchtmeyer  and  Morton, 
1967)  Indicate  that.  In  general.  It  seems  a fairly  accurate  (and  fast) 
way  of  accounting  for  shocks. 

There  are  two  places  In  our  model  where  we  expect  that  shocks  might 
form,  and  hence,  q becomes  Important.  One  Is  at  the  boundary  of  a 
condensation  where  large  pressure  gradients  form.  The  size  scale  h 
of  a knot  Is  determined  by  a balance  between  radiative  cooling  and 
conductive  heating. 


(4.22) 


in/aCt) 

The  time  scale  for  growth  of  a knot  Is  the  cooling  time. 


T - 2m 

nA  ' 


(4.23) 


In  order  to  maintain  pressure  equilibrium,  plasma  must  move  fast  enough 
to  keep  up  with  pressure  changes. 


V h/T^  fts  4 X 10^^  >/a~T^'^^  . 


(4.24) 


As  long  as  these  velocities  are  damped  quickly  (l.e.  In  a time  short 
compared  to  the  cooling  time)  then  they  remain  subsonic. 


fl  = ^ - 4 X 10®  7/^  T^'"^  , (4.25) 

4 6 

and  M < 1 for  3 X 10  < T < 5 X 10  K.  However,  since  conduction  Is  the  only 


nl/4 
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dissipative  process  present,  then  these  velocities  will  not  be  damped  because 
whenever  condensations  occur,  the  radiative  cooling  time  must  be  short 


^ 1 
; 1 


^ i 


t 

i 

ir 

r 
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compared  to  the  conductive  cooling  time.  Hence,  velocities  generated 
while  the  plasma  is  hot  will  persist  and  result  in  shocks.  At  T ~ 10^  K, 
equation  (4.24)  gives  v = 30  km  s which  is  supersonic  at  T = 3 X 10^  K. 
In  Section  5«4.2,  we  examine  a simple  analytic  model  for  the  growth  of 
condensations  and  find  that  shocks  may  indeed  form  at  the  boundary 
between  the  cold  and  hot  regions.  Our  numerical  calculations  also 

f 

support  this  result. 

Shocks  also  form  when  cool,  freely-falling  material  impacts  onto 

7 -1 

the  chromosphere.  Downward  velocities  on  the  order  of  ~ 10  cm  s 
are  observed;  these  are  well  above  the  sound  speed  at  chromospheric 
temperatures.  In  our  model,  the  artificial  viscosity  q dissipates  the 
streaming  energy  of  the  falling  material.  The  usual  dissipation  mech- 
anisms, conduction  and  viscosity  vary  as  high  powers  of  T and  are 
negligible  at  low  temperatures.  However,  the  velocities  of  the  cool 
plasma  are  obviously  damped  by  some  means  since  this  plasma  is  not 
observed  to  bounce  back  into  the  corona  after  impacting  the  chromosphere. 
We  believe  that  the  most  likely  source  of  this  dissipation  is  shocks. 

4.4  Computer  Code 

In  order  to  obtain  numerical  solutions  for  the  system  (4.3)  - (4.7), 
including  the  artificial  viscosity,  we  use  a two-step,  second  order 
difference  scheme  similar  to  the  ones  described  by  Richtmeyer  and 
Morton,  1967*  We  have  three  partial  differential  equations  to  solve 
for  the  variables  v»  S and  T.  They  can  be  written  as: 
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where  V Is  a three  component  vector  (v,  s,  T)  and  F is  a nonlinear 
function  of  V and  its  spatial  derivatives.  The  pressure  and  number 
density  are  given  as  explicit  functions  of  V by  equations  (U.6)  and 
(U.7):  hence,  they  can  be  obtained  directly.  Assuming  that  V(x,  t) 
is  known,  V(x,  t + At)  can  be  obtained  by  a Taylor  series  expansion. 

But, 
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hence 


2 


(4.28) 


(4.29) 


From  equation  (4.28)  we  note  that  F + /it  needed  only  to  first 

~2 

order  accuracy  in  At;  therefore  we  first  generate  provisional  values 

of  V K . defined  by: 

|t  + A^ 

2 


and  define 


"It+At""- 


0(At^) 


(4.31) 


Since  F is  a smooth,  well-behaved  function  of  V,  we  can  expand  the 
right-hand  side  of  equation  (4.31)  to  obtain 
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The  scheme  consists  of  first  generating  ^ from  (4.30)  and  then 

~ ~2 
calculating  ^ + ^t’  ^]t  + At  obtained  from 

'^It  + At  = ^It  ^It  + 


Equation  (4.33)  follows  immediately  from  equations  (4.29)  and  (4.32). 

The  boundary  conditions  S(x,  O),  v(x,  O)  and  T(x,  O)  are  used  to 
initiate  the  calculation  and  the  spatial  boundary  conditions  are  used 
to  compute  the  velocities  and  heat  flux  at  the  ends  of  the  space  grid. 
To  facilitate  the  computation,  we  rewrite  equation  (4.7)  as 


where 


.(X.  , 


AI(s)  = / A(s)  ds 


Three  known  functions  of  S,  g (S),  A(S)  and  Al(S),  appear  in  the 
equations  and  two  of  T,  A(T)  and  T . Rather  than  evaluate  these 
five  functions  directly  (this  requires  too  much  computer  time),  we 
first  tabulate  them  for  ~ 800  semi- logarithmically  spaced  values  of 
S and  T,  and  use  linear  interpolation  to  estimate  their  value  for  a 
particular  S or  T.  The  estimated  values  are  accurate  to  better  than 


j r 


The  fineness  of  the  x grid,  Ax,  and  the  time  increments  At  are 
determined  by  the  requirements  that  they  be  smaller  than  length  scales 
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and  time  scales  of  Interest  and  that  the  difference  scheme  remain  stable. 

The  equations  that  we  wish  to  solve  are  highly  nonlinear.  Unfortunately, 
very  little  is  known  about  the  stability  or  convergence  of  finite  difference 
schemes  for  such  systems.  The  only  accurate  method  known  to  determine 
stability  is,  of  course,  trial  and  error.  However,  we  calculate  stability 
criteria  for  several  linear  approximations  to  our  system  and  hope  that 
these  restrictions  are  sufficient  to  assure  convergence. 

If  we  neglect  the  heat  loss  terms,  the  area  factor  and  gravity, 
our  equations  reduce  to  those  of  a perfect  fluid.  Our  scheme  then 
becomes  equivalent  to  the  Lax-Wendroff  scheme  (Richtmeyer  and  Morton, 

1967)*  which  has  as  its  stability  criterion 

AS. 


f " (|v|’'i  Cj  ) 


(4.36) 


where  is  the  Eulerian  width  of  the  grid  spacing  Ax  at  the  jth  grid 


point,  and  is  the  speed  of  sound  at  this  point.  The  index  j runs 


over  all  the  intervals  in  the  x grid.  If  we  approximate  the  heat 


equation  (4.5)  by  a simple  diffusion  equation 

3T 


3t 


= K 


ds" 


(4.37) 


then  the  stability  criterion  is 


At  < Min 
j 


(As)*; 


'J 


(4.38) 


Since  the  radiative  term  is  also  Important,  we  introduce  another  restriction 


At  < Min  p /J-,,  , 

5 ^ J J 


(4.39) 


where  is  the  pressure  at  the  Jth  grid  point  and  is  the  radiative 


\ 4 
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losses  there.  Combining  these  three  criteria,  we  adjust  the  time 
increment  at  each  time  t so  that 


At  < Min  (Atj^,  Atg,  At^) 


The  spatial  grid  spacing,  Ax,  is  kept  fixed  throughout  a run  in 
order  to  assure  proper  centering  of  the  difference  equations.  From 
symmetry  arguments.  Section  4.2,  only  half  the  loop  need  be  included. 

Q 

We  use  33  points  for  the  x-grid;  hence,  for  a loop  of  height  2 X 10  cm., 
the  average  width  of  each  grid  interval  is  700  km.  The  width  of  a 
particular  interval,  i.e.  of  a particular  slab  of  plasma,  depends  on 
the  density  and  area  at  that  point;  from  equation  (4.7)  we  have: 


nA(s)  ' 

Equation  (4.4l)  implies  that  in  regions  where  the  density  is  high  or 
the  area  large,  such  as  in  a condensation  at  the  top  of  the  loop,  the 
Eulerlan  spacing  of  the  grid  is  finer.  As  mentioned  previously,  this  is 
the  main  advantage  to  using  a Lagrangean  coordinate  system. 

For  the  temperatures,  densities  and  sizes  of  typical  loops,  and 
with  32  slabs  for  the  x grid,  we  find  that  the  time  increment  At, 
as  determined  by  equation  (4.4l),  ranges  from  At  ~ .05  - .25  sec.  To 
follow  a typical  loop  from  the  10  K temperature  level  until  the  time 
when  most  of  the  plasma  has  cooled  and  fallen  onto  the  chromosphere 
requires  ~ 1 - 3 minutes  of  computer  time  on  Stanford's  IBM  370/168. 

\ large  fraction  of  this  time  is  taken  up  during  the  initial  conductive 
cooling  stage  of  the  loop  when  criterion  (4.38)  dominates.  Unfortunately, 
the  time  for  a run,  varies  as  a high  power  of  the  grid  spacing.  If 


I 


thus,  we 


_o 

equation  (4.38)  is  the  important  criterion,  then  Tj^  “ (Ax)  ; 
are  severely  limited  in  the  possible  fineness  of  the  x grid.  Details 
of  our  difference  scheme  along  with  a copy  of  the  code  are  given  in 
Appendix  B. 
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5.  COOLING  OF  FLARE  LOOPS 
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5.1  General  Evolution  of  Loops 

From  our  computer  calculations  the  following  picture  emerges  of 
the  evolution  of  a post-flare  loop.  The  evolution  can  be  divided  Into 
three  distinct  stages  characterized  by  the  behavior  of  the  temperature 
profile.  In  the  first  stage  ~ Immediately  after  flare  heating  ends  — 
the  temperature  Is  high  enough  so  that  conduction  to  the  chromosphere 
dominates  the  energy  losses.  We  expect  that  evaporation  may  also  be 
an  Important  cooling  mechanism  during  this  phase.  During  the  first 
stage,  the  plasma  Is  stable  against  the  growth  of  thermal  perturbations; 
hence,  a smooth  temperature  profile  develops.  The  scale  height  Is  so 
large  that  the  loop  Is  approximately  Isothermal,  and  the  cooling  rate 
varies  as  a high  power  of  the  temperature.  These  results  also  hold 
If  evaporation  Is  Included  (Antlochos  and  Sturrock,  19761>). 

As  the  temperature  decreases  the  effectiveness  of  conduction  drops 
rapidly  due  to  the  strong  temperature  dependence  of  the  thermal  conduc- 
tivity, equation  (3.7).  and  eventually  radiation  begins  to  dominate  the 
energy  losses.  The  loop  then  enters  the  second  stage  of  Its  development, 
during  which  radiation  Is  the  primary  cooling  mechanism.  It  Is  during 

this  phase  that  thermal  Instabilities  can  occur;  hence,  the  evolution  \ 

\ 

of  the  loop  during  the  second  stage  Is  highly  complex  involving  large  i 

temperature  gradients  and  plasma  velocities.  We  find  that  the  details  ^ 

I 

of  the  evolution.  In  particular  whether  condensations  form,  are  highly 
sensitive  to  the  Initial  configuration  of  the  loop  plasma  when  It  first 
enters  the  radiative  cooling  stage. 
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Since  the  density  in  post-flare  regions  Is  high,  radiative  losses 

are  large  and  the  temperature  throughout  a loop  quickly  drops  to  chromo- 

4 

spheric  values,  T ~ 10  K.  Once  it  reaches  this  value  the  temperature 

changes  very  little  during  the  subsequent  development  of  the  plasma 

because  the  heat  losses  at  low  temperatures  are  Insignificant.  We 

define  this  phase  of  the  evolution,  during  which  the  plasma  remains 
4 

isothermal  at  'n.  lo  K,  as  the  third  stage.  In  this  phase  large  density 
gradients  and  velocities  are  generated  due  to  pressure  forces  and  the 
effect  of  gravity. 

In  our  model,  the  end  point  of  the  evolution  Is  the  establishment 

4 

of  a gravitational  scale  height  atmosphere  at  3 X 10  K.  This  does  not 
agree  with  the  true  solar  atmosphere  since  we  have  not  Included  coronal 
heating  In  the  calculations.  For  our  large  Initial  densities,  the 
energy  loss  rate  of  the  loop  plasma  Is  much  larger  than  the  estimated 
value  of  the  coronal  heating  rate,  ^ 10°  ergs  cm  s , (Moore,  1972). 

However,  once  the  loop  plasma  cools  to  chromospheric  values  and  begins 
to  fall  onto  the  chromosphere,  then  the  density  In  the  upper  part  of 
the  loop  decreases  until  coronal  heating  once  again  becomes  Important. 

Therefore,  we  terminate  our  calculations  when  the  density  at  the  top 
of  the  loop  drops  by  more  than  an  order  of  magnitude  from  Its  Initial 
value. 

The  time  scales  for  the  three  stages  of  evolution  are  given  by: 

n 

(l)  the  conductive  cooling  time  at  high  temperatures,  T~  10'  K;  (11) 
the  radiative  cooling  time  at  T ~ 10^  K;  and  (ill)  the  free-fall  time 
along  the  loop.  For  typical  loops,  n ■ 5 X 10^*^  cm  ^ and  H » 3 X 10^  cm, 

i 

these  times  are: 


8 


(5.1) 


(11)  T,,  - 10^  8 , 


'll  2. 
n A 


(5.2) 


(til) 


10^  8 . 


(5.3) 


We  note  that  the  time  scale  for  8tage  (11)  la  shorter  than  the  llfetlass 
of  condensations  In  loop  prominences  ~ 10  minutes;  hence,  the  bright  Hor 
knots  cannot  be  due  simply  to  a temperature  difference  along  the  loop. 
Instead  they  must  represent  a density  difference.  Our  calculations  verify 
this  hypothesis.  We  find  that  due  to  the  thermal  Instability  of  the 

plasma,  large  temperature  differences  form;  but,  these  die  out  quickly 

2 

<10  s.  However,  the  temperature  gradients  give  rise  to  density 

3 

enhancements  which  persist  for  long  times  ~ 10  s,  because  the  speed  of 

sound  Is  small  at  low  temperatures.  For  example.  In  order  for  a 

^ 9 

condensation  at  10  K to  expand  to  a size  of  10  cm  at  the  sound  speed, 
the  time  required  Is  10  min. , which  Is  of  the  same  order  as  the  life- 
time of  Hot  knots. 

In  summary,  we  find  that  the  evolution  of  post-flare  loops  Is  In 
three  stages: 

(1)  conductive  cooling  (with  or  without  evaporation)  at  high  temperatures, 
~ 10^  K, 

(11)  radiative  cooling  at  lower  temperatures,  ~ 10^  K, 

(111)  the  dynamics  of  Isothermal  plasma  at  chromospheric  temperatures, 

k 

~ 3 X 10  K. 

The  stability  of  a particular  loop  Is,  of  course,  determined  by  the 
second  stage.  In  the  rest  of  this  chapter  we  describe  each  of  these 


m 


I 
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stages  In  detail,  and  determine  under  what  conditions  will  thermal 
Instability  significantly  affect  the  evolution  of  loop  plasma. 

5.2  Conductive  Cooling  Stage 
5.2.1  Conductive  Damping  Rate 

During  the  Initial  phase  of  flare  cooling  the  temperature  Is  high 
enough  so  that  the  conduction  term  dominates  In  the  dispersion  relation, 
equation  (l.l4).  Therefore,  we  expect  perturbations  to  be  damped  on  a 
time  scale. 


i •. 

I I 

• I 
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KT(f) 

since  Tjj  Is  the  time  scale  for  temperature  differences  over  lengths 


(5.4) 


X/2  to  be  significantly  reduced  by  conduction.  We  tested  this  result  by 
examining  the  behavior  of  a perturbation  of  small  wavelength  but  with 
a finite  amplitude.  The  perturbation  was  chosen  to  be  Initially 
Isobarlc  with  an  amplitude  of  9*75^  and  a wavelength  \ ® 2 X 10^  cm. 

The  loop  was  assumed  to  be  due  to  a point  dipole  source,  with  a height 
above  the  chromosphere  of  3 X 10^  ^m  and  a compression  factor  T « 10. 

For  this  choice  of  height  and  T,  the  loop  length  turns  out  to  equal 
4 wavelengths  of  the  perturbation.  The  Initial  average  temperature 
and  density  assumed  were  T * 10^  K and  n = 5 X 10^^  cm 

In  Figure  5*I>  the  evolutions  of  the  amplitudes  of  the  temperature 
and  density  variations  are  shown.  We  note  that  the  damping  time  for 
the  temperature  perturbation  Is  ~ 5 sec.,  which  Is  In  good  agreement 
with  eq.  (5.4);  however,  the  evolution  of  the  density  is  very  different 
from  that  of  the  temperature.  The  Initial  response  of  the  density 
profile  lags  behind  that  of  the  temperature  due  to  the  Inertia  of  the 
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Figure  5*1  EVOLUTION  OF  THE  AMPLITUDE  OF  A 
TEMPERATURE  AMD  DENSITY  PERTURBATION  THAT  IS 
DAMPED  BY  CONDUCTION.  The  solid  line  refers 
to  the  density  perturbation  and  the  broken 
line  refers  to  the  temperature  perturbation. 


plasma.  The  maximum  velocity  generated  by  the  resulting  pressure 


which  Is  an  order  of  magnitude  less 


than  the  sound  speed 


By  t l4  sec.  we  note  from  Figure  ^.1  that  the  amplitude  of 
the  perturbation  has  dropped  by  over  60^,  whereas,  the  thermal  energy 


of  the  loop  has  decreased  due  to  conduction  out  the  base  and  radiation 


by  less  than  1^.  The  short  damping  time  Implies  that  we  are  justified 


in  assuming  an  Isothermal  and  static  initial  state  for  the  loop  plasma 


All  small  scale  gradients  due  to  the  Initial  evaporation  of  chromospheric 


material  by  flare  heating  decay  quickly  compared  to  the  cooling  time 


It  is  possible  that  large  scale  motions,  on  the  order  of  the  length  of 


the  loop,  and  of  small  amplitude  persist  for  time  scales  on  the  order 


of  the  cooling  time.  However,  our  computer  calculations  Indicate 


that  velocities  which  are  small  compared  to  the  sound  speed  have 


negligible  effect  on  the  evolution  and  stability  of  flare  loops 


If  evaporation  cooling  Is  Important,  then  large  plasma  velocities 


may  be  driven  by  this  process.  We  can  estimate  the  magnitude  of  these 


velocities  — the  downward  heat  flux  must  be  approximately  equal  to  the 


mechanical  flux  upwards 


Taking  the  heat  flux  to  be  of  order  10  - 10^  ergs  cm  s 


and  small  loops,  respectively,  and  the  pressure 


we  obtain 
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Since  the  sound  speed  at  10  K is  ^0  km  s , the  Mach  number  is 
M K v/c  .2,  and  the  ratio  of  kinetic  to  thermal  energy  is  only 
•04.  In  addition,  we  expect  that  since  the  heat  flux  F « t , the 
Mach  number  M « T^,  if  the  pressure  is  constant,  or  M * T^,  if  density 
is  constant.  In  either  case  we  note  that  the  Mach  number  varies  as  a 
high  power  of  T;  hence,  even  if  the  evaporation  velocities  are 
comparable  to  the  sound  speed  when  the  loop  temperatures  are  high, 
these  velocities  will  be  negligible  when  the  temperature  drops  low 
enough  so  that  radiative  cooling  dominates. 


5.2.2  Static  and  Evaporative  Cooling  Rates 


The  first  stage  cooling  of  flares  has  been  investigated  by  several 
authors;  Culhane  et  al.  (1970),  Strauss  and  Papaglannis  (1971),  and 
Zaumen  and  Acton  (1974).  In  their  models  they  assumed  that  a loop 
has  uniform  cross-section;  however,  Antlochos  and  Sturrock  (1976a) 
point  out  that  for  magnetic  fields  typical  of  flare  loops  (Rust  and 
Bar,  1973 )»  the  variation  of  area  along  the  loop  has  a significant 
effect  on  the  cooling.  Antiochos  and  Sturrock  have  calculated  the 
evolution  of  the  temperature  and  density  profiles  both  for  static  and 
evaporative  (1976b)  models;  they  assume  that  gravitational  and  radiative 
effects  are  negligible  and  that  all  velocities  are  small  so  that 
pressure  gradients  are  unimportant.  With  these  assumptions  the  MHD 
equations  simplify  sufficiently  so  that  analytic  solutions  may  be 
obtained. 

For  completeness  we  state  the  important  results  of  the  calculations 
by  Antiochos  and  Sturrock  below.  In  the  case  of  static  cooling  of  a 
dipole  loop,  the  temperature  profile  is  given  by 


N 


where  R is  the  radius  of  the  loop.  Figure  (3«1).  The  parameter  is 
directly  related  to  the  compression  factor  F — the  dependence  of  t on 
F is  given  in  Figure  2 of  Antiochos  and  Sturrock-  (1976a).  For  the 
case  of  evaporative  cooling,  the  results  are  similar  except  that  the 
temperature  varies  as  (1  + t/r)  and  at  low  temperatures  the 
heat  flux  varies  linearly  with  T rather  than  being  constant  as  in  the 
static  case. 

Q 

We  ran  our  code  for  a dipole  loop  with  F • 10  and  R = 2.5  X 10  cm, 
in  order  to  compare  our  model  with  the  results  above.  The  initial 
temperature  and  density  were  chosen  as  1.05  X 10^  K and  2.5  X 10^^  cm 


respectively.  With  these  values  the  cooling  time  given  by  equation 
(5.8)  is  200  sec.  In  Figure  5*2  we  plot  ( ^(o))  “ ^ versus  t, 

where  T is  the  average  temperature  of  the  loop  in  our  computer  run. 


We  note  that  the  graph  is  linear  in  agreement  with  equation  (5*7)  and 


chat  the  cooling  time  is  170  sec.  which  also  agrees  with  theory.  The 


loop  size  and  plasma  density  were  deliberately  chosen  to  be  small  so 
chat  radiative  cooling  would  be  negligible.  We  find  that  Che  plasma 
loses  32')^  of  its  thermal  energy  in  270  sec.,  with  88^  of  Che  loss  due 
Co  conduction  and  the  rest  via  radiation. 
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These  results  indicate  that  in  the  limit  of  small  radiative  losses 
and  plasma  velocities  our  numerical  code  does  reproduce  the  analytic 
calculations.  Of  course,  this  cannot  be  used  as  an  argument  to  prefer 
static  over  evaporative  cooling  since  our  model  assumes  no  evaporation. 

In  fact,  we  expect  that  evaporation  is  probably  the  dominant  cooling 
process  at  high  temperatures  because  the  heat  flux  predicted  by  both 
our  code  and  the  analytic  model  is  high.  The  initial  heat  flux  is 

Q -2-1 

2 X 10^  ergs  cm  s , which  is  much  too  large  to  be  absorbed  by  the 
chromosphere.  We  believe  that  the  initial  cooling  is  accurately 
described  by  the  model  of  Antiochos  and  Sturrock  (1976b)  and  that  the 
later  stages,  when  radiation  is  important,  are  accurately  described  by 
our  numerical  code. 

From  the  arguments  in  the  previous  two  sections,  we  reach  the 
following  conclusions.  The  evolution  of  the  loop  density  and,  hence, 
the  stability  of  the  plasma,  are  not  uniquely  determined  by  the  Initial 
temperature  and  density  since  evaporation  may  occur.  Even  if  evaporation 
is  present,  however,  the  first  stage  cooling  process  will  bring  the 
loop  into  an  approximately  isothermal >and  static  state.  Given  the 
density  at  lower  temperatures,  when  radiation  losses  dominate, 
the  stability  may  be  determined.  In  the  next  section  we  discuss 
necessary  criteria  on  the  density  and  temperature  profiles  for  the 
formation  of  condensations. 

5.3  Instability  Criteria 

Using  our  model  of  flare  loops,  we  have  determined  criteria  under 
which  thermal  instabilities  in  a coronal  plasma  can  grow.  The  term 
instability  is  not  strictly  applicable  to  our  model  since  the  plasma 
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is  never  in  an  equilibrium  state  but  is  continuously  evolving;  however, 
we  use  it  to  denote  situations  in  which  the  cooling  is  highly  non- 
uniform.  By  the  term  "unstable  loop",  we  mean  that  a condensation 
observable  in  Ho  and  of  small  size  scale  forms.  Specifically,  we 
require  that  the  condensation  be  at  least  an  order  of  magnitude  brighter 
in  Ho  than  the  rest  of  the  loop  and  have  a lifetime  of  the  order  of 
minutes.  These  condensations  correspond  to  the  bright  knots  in  HO 
that  are  often  observed  to  form  at  the  tops  of  loop  prominences 
(Section  2,2.4). 

We  find  that  the  following  conditions  are  necessary  for  a loop 
to  be  unstable  (according  to  the  definition  above): 

(1)  Since  both  the  growth  time  of  thermal  perturbations  and  the  lifetime 
of  a loop  are  of  the  order  of  the  radiative  cooling  time,  only  pertur- 
bations of  finite  amplitude  can  lead  to  observable  effects.  However, 

we  find  that  the  amplitude  required  is  quite  small,  of  order 

(2)  From  linear  theory  (Section  1.2),  we  require  that  the  growth  rate 

due  to  radiative  cooling  dominates  the  damping  rate  due  to  conduction. 
Hence  from  equations  (1.19)  (1.20 ) we  require. 


nX  4 X lO"^  cm"^  , (5.9) 

where  X is  the  wavelength  of  the  perturbation  and  A is  the  radiative 
loss  function. 

(3)  In  order  for  a density  perturbation  to  grow,  the  cooling  time  must 
be  long  compared  to  the  speed  of  sound  travel  time  across  the  pertur- 
bation; hence  from  equation  (1.23), 

nX  « 4 X lO"^^  cm"^  . (5.10) 
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(k)  Additionally,  we  find  that  for  temperatures  and  size  scales  relevant 
to  loop  prominences,  the  number  density  in  a loop  must  be  greater  Chan 
approximately  10^^  cm  For  densities  less  than  this  value  the  radiative 
cooling  time  is  so  long  Chat  gravitational  effects  and  conduction  to 
the  chromosphere  dominate  the  cooling  process,  stopping  the  growth  of 
condensations.  Also,  for  low  densities  coronal  heating  becomes  important 
and,  hence,  our  model  is  no  longer  accurate. 

The  conditions  enumerated  above  agree  well  with  observations  that 
loop  prominences  occur  only  after  flares  when  the  coronal  density  is 
enhanced  and  only  on  large  loops  of  size  scales  ^ 10^  cm.  In  Figure 
(5*3)  we  graph  the  relations  (2)  and  (3).  Inequality  (2)  is  valid 
in  the  region  to  the  left  of  the  broken  line  labeled  (ii)  and  inequality 
(3)  is  valid  to  the  right  of  the  line  labeled  (iii).  In  the  region 
between  line  (ii)  and  (iii)  isobarlc  perturbations  will  grow.  If 
conditions  (1)  and  (4)  are  satisfied  as  well,  then  observable  conden- 
sations will  result. 

5.4  Radiative  Cooling  Stage 

5.4.1  Comparison  with  Linear  Theory 

We  define  the  second  stage  of  the  cooling  process  to  be  that  phase 
during  which  radiation  is  the  dominant  energy  loss  mechanism.  At  the 
onset  of  instability,  we  expect  that  the  linear  theory  may  approximately 
describe  the  evolution.  In  order  to  compare  our  nonlinear  model  with 
this  theory,  we  followed  with  our  code  the  evolution  of  an  initial 
perturbation  in  the  temperature  and  density  profile  of  a loop.  The  loop 
was  chosen  in  each  case  to  be  that  due  to  a dipole  field  with  a compression 
factor,  r » 10.  The  perturbation  was  chosen  in  each  case  to  be  initially 
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isobaric,  with  an  amplitude  of  10^  and  a wavelength  equal  to  the  length 
of  the  loop,  l.e.  we  assumed  as  Initial  conditions; 

T(x,0)  = T(0,0)  jl  - .1  cos^-^jj  (5.11) 

and 

p(x,0)  = p(0,0)  , (5.12) 

where  ± Xj„  are  the  Lagrangean  positions  of  the  base  points  of  our 
mode  1 loop . 

From  conditions  (5«9)  (5.10),  the  important  parameters  that 

characteri?e  the  stability  of  a loop  are  the  temperature  and  the  total 
number  of  electrons,  nX,  where  X is  the  wavelength  of  the  perturbation. 

In  each  of  the  runs  the  values  of  T(0,0)  and  nX  were  chosen  so  that 
condition  (5.9)  was  definitely  not  satisfied  and  the  perturbation  was 
initially  damped  by  conduction.  However,  as  the  plasma  cooled  and  T 
decreased,  (5*9)  eventually  became  satisfied  and  at  a critical  temperature 
T^,  the  amplitude  began  to  Increase.  Figure  5.3  shows  the  value  of 
as  a function  of  nX.  There  is  some  ambiguity  in  the  exact  value  of 
since  in  no  case  did  the  temperature  profile  retain  the  simple  form 
given  by  equation  (5.11).  Instead,  conductive  losses  to  the  chromosphere 
and  small  but  non-negligible  velocities  altered  the  profile.  Therefore, 
we  assumed  that  the  temperature  and  density  variations  were  due  to  two 
parts , 


^ + (^\  I 

T~\Tja  \T/i  / 

An  /An\  1 

n^ln/a  \n/i  / 
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where 


Hi  - H\ 


Indicate  variations  due  to  adiabatic  motions  of  the 


plasma  and  indicate  the  amplitude  of  the  isobaric  pertur- 


bation. Since  for  adiabatic  motions 


l^\  2 /^\ 

\ Tfa"  3 \ n/a  ’ 


(5.14) 


the  perturbation  amplitude  is  ^iven  by 


{^\  . I^\  . i ^ 2 ^ 

\T/i“"\n/i*5  t"5  n * 


(5.15) 


Hence,  we  define  as  the  average  temperature  In  the  loop  at  which  the  value 


as  given  by  equation  (5.15)»  began  to  increase. 


In  general  the  results  of  our  computer  calculations  agree  well  with 
theory.  We  find  that  is  approximately  Independent  of  n or  X as  long 

as  the  product  nX  is  constant.  From  Figure  ^.3  we  note  that  the  effect 
of  the  nonlinearities  in  our  model  is  to  Increase  the  stability  of  the 
loop  especially  at  high  temperatures.  This  Is  to  be  expected  since  for 
the  initial  perturbation  with  the  longest  possible  wavelength,  and, 
hence,  the  fastest  growth  rate  (equation  ^.11),  the  temperature  maximum 
occurs  at  the  base  of  the  loop.  Conduction  cooling  to  the  chromosphere 
decreases  the  temperature  near  the  base,  thereby  decreasing  both  the 
amplitude  and  the  length  scale  of  the  temperature  variations.  We  find  that 
loops  with  a large  compression  factor  F can  support  a perturbation  with  a 
relatively  larger  wavelength  because  heat  losses  to  the  chromosphere  are 
reduced;  hence,  they  are  more  susceptible  to  Instability,  l.e.  In- 
creases with  r.  Once  instability  sets  in,  however,  the  evolution  of  the 
plasma  is  insensitive  to  F. 

In  all  cases  we  find  that  pressure  balance  is  maintained  for  only 
a short,  initial  part  of  the  cooling.  As  can  be  seen  in  Figure  ^.3,  the 


temperature  occurs  near  the  line  labelled  at  which  condition  (^.10) 
breaks  down  and  the  perturbation  can  no  longer  remain  Isobaric.  There- 
fore, significant  velocities  soon  develop.  In  the  next  section  we 
Investigate  the  growth  of  these  velocities  using  an  analytic  model  of 
a growing  condensation. 

5.4.2  Analytic  Model  of  a Condensation 

We  have  constructed  a simple  analytic  model  to  follow  the  growth 
of  a condensation  into  the  nonlinear  regime.  The  model  is  valid  during 
the  initial  growth  when  pressure  equilibrium  is  maintained,  i.e. 

^p  « 0 . (5 .16) 

Equation  (5.16)  will  be  valid  as  long  as  the  speed  of  sound  and  the 
gravitational  scale  height  are  large  compared  to  relevant  velocities 
and  size  scales.  This  model  can  be  used  to  calculate  the  temperature 
at  which  velocities  resulting  from  the  difference  in  cooling  rates 
between  the  cold  and  the  hot  plasma  become  appreciable  and,  hence, 
pressure  equilibrium  breaks  down.  In  addition,  we  can  obtain  estimates 
of  the  magnitude  of  this  velocity  and  determine  if  shocks  form.  Of 
course,  these  results  can  also  be  obtained  using  our  full  computer 
code;  however,  we  find  that  the  analytic  model  is  useful  both  to 
check  the  computer  calculation  and  to  clarify  the  important  features 
of  the  nonlinear  instability. 

We  assume  that  the  temperatures  are  low  enough  so  that  heat 
conduction  is  negligible  compared  to  radiative  losses.  Therefore  the 
hot  and  the  cold  regions.  Figure  5*4,  may  be  taken  to  be  separated  by 
a boundary  of  negligible  width.  Within  each  region  the  temperature 
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and  density  are  assumed  to  be  approximately  uniform  and  given  by, 
(l^(t),  and  |Tj^(t),  where  the  subscripts  "c"  and  "h" 

refer  to  the  cold  and  hot  plasma  respectively.  Since  the  two  regions 
cool  radlatlvely  at  different  rates,  the  position  of  the  boundary  must 
change  In  order  to  maintain  pressure  balance,  equation  (^.l4);  however, 
since  conduction  Is  neglected  no  material  moves  across  this  boundary. 

From  Figure  5»4  we  note  that. 


L (t)+i  (t)=L/2»  const.  , 
c ' ' n 


thus  the  volume  of  the  condensation  Is  given  by: 


v,(t)  - sy^‘ 


A(s)  ds  , 


and  of  the  hot  plasma: 


Vh(t)  . - v^(t)  , 


where  V^,  the  total  volume  of  the  loop.  Is  constant.  Differentiating 
equation  (^.18)  we  obtain  the  rate  of  change  of  the  volumes: 


^ Vjt)  . - V^(t)  - 2A(i^)v  . 


dt  h' 


where  v is  the  velocity  of  the  boundary, 


Since  the  temperature  and  density  are  uniform  over  each  region,  we  can 


Integrate  the  mass  equation  over  the  volume  of  the  condensation: 
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X It  V . V d3v  - 0 . 


Using  the  divergence  theorem  on  the  second  term  of  equation  (5«22)  and 
substituting  equation  (^.20)  yields, 


therefore , 


"cCt)  V^(t)  - const.  (5.2l 

Equation  (^,2h)  is  equivalent  to  the  assumption  that  there  Is  no  mass 
transfer  between  t^ie  hot  and  cold  regions.  This  can  only  be  true  If 
the  cold  plasma  cools  faster  than  the  hot,  l.e.  the  temperature  must 
be  In  the  unstable  regime  of  the  cooling  curve.  Once  the  temperature 
enters  the  stable  regime,  the  profile  must  eventually  become  uniform 
throughout  the  loop. 

The  heat  equation  can  be  Integrated  In  a similar  manner  as  the  mass 
equation;  yielding 


^ + ^p  — V 
2 c dt  2 P dt  c 


•vo’'r''c  ■ 


where  we  have  assumed  a radiative  loss  function  of  the  form. 


f(n,T)  . AqhV® 


Defining, 


\i^)  • vjt)/v^ 
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\(t)  = Vj^(t)/V^  , 


(5.27b) 


we  obtain  the  following  set  of  equations; 


X^(t)  + \{t)  = 1 

= ”c,h  = 

o ujfc  Cjh  ■ - A^n 


2 c.h  dt  2 dt 


„2  _,-0f  , 

0 c,h  c,h  c,h 


P “ 2kn^.h\.h  j 

The  seven  Independent  equations  In  the  system  (5.28)  are  enough  to 
uniquely  determine  the  variables  (n^  \^),  (n^^  \ P 

given  Initial  conditions. 

The  quantity  that  is  of  interest  Is 


5(t)  = Tjt)/Tj^(t) 


the  ratio  of  cold  to  hot  temperature.  (The  amplitude  of  the  perturbation 
Is  equal  to  (1  - S)/(l  -k  S).)  After  some  manipulation,  the  two  heat 
equations  In  (3.28)  can  be  decoupled  and  integrated  to  yield. 


F(5)  = F(^)  - P(S^)t/T  . 


where 


1 _S(QH-1)  2C»-1 

F(?)  *y  1 - (^+  S)  ^ d§  , 

r -|  5(eH-l)  2a- 1 

- lb  - C * 
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''o  "c<0) 


(5.33) 


n - H„/N, 


n 

' ( 
f 

i 


There  are  two  adjustable  parameters  In  our  solution,  T|,  the  ratio  of  the 
total  number  of  cold  particles  to  hot  and,  the  Initial  temperature 
ratio  S(0).  We  are  Interested  In  Initial  perturbations  of  small 
amplitude  and  large  wavelengths;  hence. 


• 9 ^ ? <1  and  r|w  1 . 

o 


For  example,  an  Initial  perturbation  with  an  amplitude  of  and  a 
wavelength  equal  to  the  length  of  the  loop  corresponds  to  ^ > .9  and 
Ti-  1. 

In  deriving  equations  (^.30)  and  (^.31 )>  we  assumed  that  the 
exponent,  or,  of  T In  the  radiative  loss  function  (equation  (3«26)),  is 
a constant.  From  Cox  and  Tucker's  (I969)  results.  Figure  1.1,  we 
note  that  or  Is  actually  a function  of  T;  however,  we  can  approximate 
the  cooling  function  by  segments  with  constant  01  as  In  Moore  and  Fung 
(1972)  and  Goldsmith  (1971).  Since  the  cooling  function  Is  probably 
only  accurate  to  at  best  a factor  of  two,  the  approximation  of  constant 
or  will  not  affect  the  general  results  of  our  model. 

The  time  dependence  of  $ can  be  determined  from  equation  (5>30)  and 
(3*31. )i  however,  they  break  down  for  or  ■ «2  as  can  be  seen  from  the  form 
of  (5*31  )•  This  Is  to  be  expected  since  we  assumed  that  the  loop  Is 
unstable  to  the  growth  of  the  condensation,  but  for  a cooling  function 
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of  the  form  given  by  equation  (5.26),  Field's  (I965)  criterion. 


n ax  ^ , 
dT  ■ T 3n 


(5.36) 


is  not  satisfied  for  Ot  < -2.  Thus,  our  solution  Is  only  valid  for 
a > -2.  For  certain  a equation  (5.31)  can  be  integrated  directly; 


F(?) 


(5.37) 


From  the  Cox  and  Tucker  curve  we  find  that  1/2  ^ 2 in  the  unstable 

5 7 

temperature  regime,  10  ^ T ^ 10  . 

In  Figure  (5.5)  we  plot  g as  a function  of  t/r  for  each  of  the 
three  values  of  a in  equation  (5.37)  and  for  Tj*  1,  !^  = .9.  We  note 
from  Figure  (5.5)  that  for  g < lo”^,  g is  a very  rapidly  decreasing 
function  of  t.  Since  the  plasma  tries  to  maintain  pressure  balance 
between  the  hot  and  cold  regions,  we  will  find  that  the  sharp  decrease 
in  g implies  equally  rapid  changes  in  the  density  and  hence,  supersonic 
velocities  result.  The  results  in  Figure  (5.5)  afc  not  sensitive  to 
the  value  of  T)  or  Defining  t^  as  the  value  of  t when  g * 0,  we 

find  from  equation  (5.30),  (5.32)  and  (5.37)  that  for  a ■ 1/2,  t^  * 1, 
independent  of  Ti  or  g^.  For  » » 2,  t^  -•  const,  as  T) -•  ® or  0,  or  as 
Fq  -•  1.  For  a ■ -1,  tj  diverges  for  r| -•  0 or  g^  -•  1,  but  only  logarith- 
mically. (Also,  we  are  mainly  interested  only  in  positive  values  of  a.) 

To  obtain  the  time  dependence  of  the  other  variables , we  express 
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them  in  terns  of 


\ * 5/(n+  5 ) . 


["*  ^)  (i  - CV 


T /T  - S. 

*=  Pq  ^0 


n /a  . = \ ^/\ 
c cO  cO  c 


The  equations  for  and  n^  are  identical  to  (5«^0)  and  (5*^1)  except 
that  "c"  is  replaced  by  "h". 

Our  model  breaks  down  when  the  velocities  generated  by  the  unequal 
cooling  rates  become  comparable  to  the  speed  of  sound.  We  calculate  v, 
the  velocity  of  the  boundary  between  the  two  regions.  Using  equations 
(5.20),  (5.27)  and  (5.38)  yields, 

L <A>  n d?  , 


where  we  have  defined. 


2 \ 

<A>  = Jy  A(s)  ds  = ^ . 


Defining, 


th  2t  ’ 


and  using  equations  (5.30)  - (5.32)  we  obtain. 


I 
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8 QH-11  p ^ 

2 o>  [1  - C 

5 A(A  ) 5 COH-1 ) 20H^ 

[,  . j«2]3(0H.2)  J,  3 


(S.ltS) 


If  the  loop  is  of  coasCant  cross-section  then,  of  course. 


A(A^) 


In  general,  we  expect  this  factor  to  be  less  than  unity  since  A increases 
as  (and  hence  5)  decreases.  For  example,  a dipole  loop  with  a 
compression  factor  F = 30  and  a volume  for  the  condensed  region,  ■ l/lO, 
implies  -j  * .46.  Therefore,  the  variable  area  factor  acts  to  decrease 
the  velocity,  improving  the  validity  of  our  model.  A lower  limit  on 
the  size  of  the  area  ratio  is 


A(Jl  ) A(s  * O) 


(5.46) 


however,  this  ratio  is  usually  much  larger  than  F . In  the  subsequent 


discussion  we  assume  that  -r-p-: — r fts  1. 

A(^j.) 


We  have  plotted  in  Figure  5.6  l~^l  as  a function  of  5 for  = 1 


and  = .9,  and  for  several  values  of  Of.  Assuming  typical  loop  parameters, 

6 10  Q c 

T ~ ^ X 10  K,  n~5>^10  cm^,  and  L ~ 10^  cm  implies  that  10; 

'^th 

hence,  our  model  will  begin  to  break  down  at  approximately  this  value  for 

V z'  V 

. From  Figure  5.6  we  note  that  = 10  corresponds  to  a value  of 

^th  ^th 

^ ranging  from  .4  - .2  for  .5  ^ or  £ 2.  These  values  of  $ correspond  to 
a perturbation  amplitude  of  43^  “ 6j%.  The  temperatures,  T^  and  T^^,  can 
be  calculated  for  a particular  5 using  equations  (5.38)  - (5.4o).  In 
Figure  5.7,  ^h^^hO  plotted  as  functions  of  | for  tl  » 1 and 

^ * .9.  note  that  for  the  values  of  ^ at  which  our  model  begins  to 
break  down,  the  condensation  temperature  has  decreased  by  approximately 
a factor  of  6 from  its  original  value  while  the  hot  temperature  has 
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dropped  by  a factor  of  2.  Since  the  plasma  must  cool  over  several  decades 
in  temperature  (from  10 ' K to  10  K),  we  expect  that  large  velocities 

(and  probably  shocks)  will  be  present  during  most  of  the  cooling. 

It  is  instructive  to  examine  the  dependence  of  these  results  on 
the  initial  amplitude  From  equations  (5*38)  - (5*^0)  we  obtain. 


’+2  v3(QH-2) 


(5.47) 


Letting  6 = 1 - and  « = 1 - 5.  equations  (5*45)  and  (5.4?)  imply  that 
as  6 -•  1, 


and 


j-3(0H-2)/5  5 


2 <A>  n 

5 (^1)2 


(oH-2) 6. 


(5.48) 


Therefore,  if  the  initial  perturbation  amplitude  6 is  small  enough,  signif- 
icant temperature  differences  or  velocities  do  not  form  even  if  the  plasma 
cools  considerably,  l.e.  T^/T^q  is  very  small.  This  agrees  with  our 
instability  criterion  (l)  that  only  perturbations  of  finite  amplitude 
can  grow  to  observable  condensations  (Section  5.3). 

In  summary,  we  find  that  the  analytic  model  predicts  large  velocities 
at  the  boundary  between  a condensation  and  the  hot  plasma  in  a loop.  The 
growth  rates  of  the  temperature  differences  and  the  velocity  are  given  by 
•quetlona  (5.30)  and  (5.45),  and  are  plotted  in  Figures  5.5  and  5.6.  For 
• irlr  SMll  perturbation  amplitude,  ~ 5^,  the  velocity  becomes  comparable 
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As  discussed  in  the  previous  section,  large  velocities  soon  build 
up  due  to  the  non-uniform  cooling.  With  our  code  we  can  follow  the 
details  of  the  temperature,  density,  and  velocity  profiles.  In  this 
section  we  Investigate  the  evolution  of  an  unstable  loop  throughout 
Its  radiative  cooling  stage.  Although  we  only  discuss  the  results  for 
a particular  loop  In  this  section,  the  general  features  of  the  evolu- 
tion turn  out  to  be  the  same  for  all  unstable  loops.  The  only  sensitive 
adjustable  parameter  Is  the  amplitude  of  the  Initial  perturbation;  we 
describe  Its  effect  on  the  evolution  In  the  following  section. 

The  temperature  and  density  profiles  for  a typical  unstable  loop 
at  different  times  of  Its  evolution  are  shown  In  Figures  5.8  and  5«9. 

The  Initial  profiles  were  assumed  to  be  of  the  form; 


T(0,x)  = T(0,0) 


.05  cos 


(el) 


(5.49a) 


n(0,x)  - n(0,0)  | 1 + .05  cos 


ce)) 


(5.49b) 


where  x Is  the  Lagrangean  coordinate  and  ± x are  the  positions  of 

in 

the  foot-points  of  our  mcdel.  Equations  (5.49a)  and  (5.49b)  Imply  that 
the  perturbation  Is  approximately  isobarlc,  centered  at  the  top  of  the 
loop  (x  ■ 0),  and  of  wavelength  equal  to  the  loop  length.  Since  the 


conduction  damping  rate  varies  Inversely  as  the  square  of  the  wavelength, 
equation  expect  perturbations  of  the  maximum  possible  wave- 


I 

* 

¥ 

s 
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length  to  be  most  susceptible  to  thermal  Instability.  The  Initial  values 
were:  T(0,0)  = 4 X 10^  K and  n(0,0)  * 5 X 10^^  cm  The  loop  was 

taken  to  be  that  due  to  a dipole  field  with  a compression  factor  of 
r = 10  and  a total  length  of  80,000  km,  which  Implies  a maximum  height 
above  the  chromosphere  of  30,000  km.  These  parameters  were  chosen  so 
that  T(0,0)  T^,  the  Initial  temperature  f o * the  onset  of  Instability 

(Section  5*4.1). 

From  Figure  5*8  we  note  that  by  t = 480  sec  the  amplitude  of  the 
temperature  variations  has  grown  significantly;  the  temperature  of  the 
hot  plasma  has  dropped  to  1/3  of  Its  original  value,  and  the  ratio  of 
cold  to  hot  temperatures,  § = .37*  However,  the  ratio  of  minimum  to 
maximum  density  (Figure  5*9)  at  this  time  Is  only  .73;  hence,  the 
temperature  perturbation  has  begun  to  run  away  from  the  density  variations. 

The  plasma  velocities  are  still  small,  the  maximum  velocity  Is  - 1.5  X 

6*1  “1 
10°  cm  s (directed  up  the  loop)  which  Implies  a Mach  number,  M » 10  . 

We  find  that  the  condensation  now  enters  a very  rapid  cooling  stage, 

so  that  I decreases  drastically.  This  behavior  Is  In  agreement  with  our 

analytical  calculations  In  the  previous  section.  By  t = 490  sec  the  top 

4 

5,000  km  of  the  loop  have  dropped  to  chromospheric  temperature,  3 X 10  K. 
At  this  time  the  loop  top  becomes  visible  In  VO.  Comparing  the  temper- 
ature profile  at  t > 480  and  t 490  sec,  we  see  that  It  requires 
< 10  sec  for  the  condensation  to  cool  from  ~ 5 X 10^  K to  Its  minimum 
temperature.  (We  fix  the  temperature  so  that  it  cannot  drop  below 
3 X 10^  K — Section  3*3)*  The  rapid  cooling  is  due  to  the  strong 
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Flgur«  5.9  DENSITY  PROFILE  OF  AN  "UNSTABLE"  LOOP 
AT  VARIOUS  TIMES  IN  ITS  EVOLUTKW.  (For  clarity 
tha  prof 11a  at  t > U80  sac.  Is  Indicated  by  a 
broken  llna.) 


radiative  losses  at  low  temperatures;  the  Cox  and  Tucker  loss  function 

peaks  at  about  2 X 10^  K.  For  T - 5 X 10^  K and  n - 5 X 10^*^  cm"^,  the 

theoretical  cooling  time  t = ^ ^ fa  1 sec. 

n A 

On  the  other  hand,  the  density  does  not  change  greatly  between  the 

times  t = 480  to  t = 490  sec.  Figure  5.9.  Hence,  large  pressure  gradients 

result.  This  Is  to  be  expected  since  the  velocities  are  too  small  to 

move  a large  mass  of  plasma  In  only  10  sec.  The  analytic  model  predicts 

-2 

that  for  these  small  values  of  ?,  < 10  , the  velocities  should  be 

highly  supersonic;  however,  that  model  neglects  Inertial  effects  while 

our  computer  calculations  indicate  that  they  are  important.  The  maximum 

velocity  at  t = 490  sec  is  - 4.2  X 10  cm/ sec,  which  is  only  slightly 

supersonic,  but  it  is  a factor  of  3 increase  over  the  value  at  t * 480  sec. 

Therefore,  while  the  evolution  of  the  loop  is  far  from  isobaric,  it  is 

also  not  isochoric;  instead,  we  find  that  the  dynamics  of  the  plasma 

is  critical  in  determining  whether  a long-lived  condensation  can  result. 

4 

As  the  plasma  continues  to  cool,  the  region  at  3 X 10  K continues 

to  expand  until  the  complete  loop  reaches  this  temperature.  This  occurs 

at  t > 330  sec.  We  note  that  it  takes  about  90  sec  from  the  time  that 

4 

the  top  of  the  loop  first  drops  to  3 X 10  K until  all  the  plasma  cools 
to  this  value.  If  the  loop  were  of  constant  density,  then  in  HQf  the 
evolution  of  a condensation  would  appear  to  follow  the  temperature 
evolution.  Since  the  length  of  the  loop  is  80,000  km,  this  would 
Imply  an  average  expansion  velocity  of  the  condensation  at  4^0  km/ sec, 
which  Is  larger  than  observed  velocities  (Section  2.2.4).  However,  this 
Is  not  the  case  because  the  density  profile  does  not  remain  uniform. 
Instead,  large  variations  develop  as  a result  of  the  non-uniform  cooling. 
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We  find  that  the  maximum  upward  velocity  during  this  phase  of  the 

8 

cooling  occurs  at  t - 520  sec  at  a point  9 X 10  cm  from  the  top  and  is 

equal  to  - 7«1  X 10^  cm/sec.  From  Figure  5*8,  we  can  see  that  the 

position  of  this  velocity  is  at  the  boundary  of  the  condensation.  The 

speed  of  sound  in  the  cold  region  is  3 X 10^  cm/sec;  hence,  a weak 

shock  occurs  at  the  boundary,  M 2.  The  boundary  density  is  approximately 

U X 10^*^  cm  therefore,  the  mechanical  energy  flux  into  the  condensation 
6 "2  "1 

is  of  order  lo'^  ergs  cm  sec  . However,  the  radiative  loss  rate  for 

2 -1  “3  “1 

the  densities  in  the  condensation  Is  so  high,  n A ^ 10  ergs  cm  sec  , 

that  the  mechanical  energy  flux  is  easily  dissipated  without  significant 

effect  on  the  temperature  of  the  condensation. 

In  Figure  5»9»  we  show  the  density  profile  at  t = 570  sec,  when  the 

4 

temperature  profile  Is  uniform  at  3 X 10  K.  As  a result  of  the  dynamics, 

a large  density  gradient  has  developed  at  the  top  of  the  loop.  The 

10  -3 

density  over  most  of  the  loop  volume  Is  approximately  3 X 10  cm  ; 

g 

however,  a small,  < 10°  cm,  clump  of  plasma  has  formed  at  the  top  with 

densities  an  order  of  magnitude  higher.  The  clump  contains  2'^%  of  the 

2 

total  mass  of  loop  plasma.  Since  the  radiative  losses  vary  as  n , the 
compression  at  the  top  would  appear  to  be  as  much  as  2 orders  of  magnitude 
brighter  in  HOT  than  the  sides  of  the  loop  (assuming  the  loop  is  observed 
at  the  solar  limb). 

We  believe  that  this  compression  corresponds  to  the  bright  knots  In 
Ho  observed  at  the  beginning  of  a loop  prominence  event.  Section  2.2.4. 

The  time  scale  for  temperature  variations,  < 90  sec.  Is  too  short  to 
account  for  the  long-lived  knots,  ^ 10  min.  But,  the  time  scale  for  the 
density  variations  does  turn  out  to  be  of  the  order  of  minutes.  Section 
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Hence,  we  find  that  even  though  the  HOT  bright  knots  are  a result 
of  the  thermal  Instability  of  coronal  plasma,  they  are  not  a direct 
manifestation  of  this  process.  Instead,  the  knots  are  mainly  due  to 
the  compression  of  loop  plasma  by  velocities  generated  during  the 
non-uniform  cooling.  As  a result,  their  size  scale  is  not  set  by  the 
wavelength  of  the  initial  perturbation,  which  is  of  the  order  of  the 
length  of  the  loop,  but  is  determined  by  the  details  of  the  plasma 
dynamics . 

The  results  described  in  this  section  are  relatively  insensitive 
to  the  particular  parameters  of  a loop.  The  main  criterion  that  must 
be  satisfied  by  a specific  perturbation  in  order  for  it  to  develop 
into  an  observable  knot  is  that  its  amplitude  at  the  onset  of  Instability 
be  ^ The  importance  of  the  plasma  dynamics  follows  from  the  fact 

that  the  cooling  time  is  much  less  than  the  speed  of  sound  travel  time 
across  the  loop  at  low  temperatures.  For  the  plasma  in  a condensation 
to  be  compressed  significantly,  the  velocity  generated  by  the  cooling 
must  be  supersonic.  If  it  is  subsonic  then,  of  course,  the  plasma 
streaming  energy  is  less  than  the  thermal  energy;  hence  large  pressure 
(density)  enhancements  will  not  form. 

We  can  estimate  the  dependence  of  the  velocities  on  the  density  and 

wavelength  of  the  Initial  perturbation.  The  velocity  generated  by 

2 

pressure  gradients  of  order  p/X,  acting  on  a time  scale,  p/n  A,  is: 

2 

V ^ -^6—  • (5.50) 

m^n^AX 

Letting  T^  be  the  temperature  at  instability  onset,  equations  (^.^0) 

and  (5.9)  yield 
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10  K,  we  can  approximate 


hence 


Equation  (^.^2)  implies  that  the  velocity  generated  by  the  non-uniform 


cooling  and,  therefore,  the  formation  of  a condensation  is  only  weakly 


dependent  on  n or  X.  Our  computer  calculations  confirm  this  result 


Effect  of  Perturbation  Amplitude 


In  agreement  with  the  analytic  model,  we  find  that  the  evolution 


of  the  loop  plasma  is  strongly  dependent  on  the  amplitude  of  the 


thermal  perturbation  at  the  onset  of  instability,  T . If  the  amplitude 


is  too  small,  then  large  temperature  differences  do  not  last  for  a 


long  enough  time  to  produce  significant  density  variations 


This  is  illustrated  by  Figures  3«10  and  ^.11  in  which  are  plotted 


the  temperature  and  density  profiles  of  a loop  with  the  same  parameters 


as  the  ones  used  in  the  previous  section,  except  that  the  initial 
amplitude  has  been  halved  to  2.^'^.  Large  temperature  differences  occur 


K while 


the  sides  are  at  10  K.  However,  the  temperature  variations  last  for 


less  than  40  sec.  By  t > ^30  sec,  the  temperature  throughout  the  loop 


Comparing  Figures  5*8  and  5>X0,  we  note  that  similar  temperature 
profiles  occur  in  the  loop  with  the  3%  amplitude  and  in  that  with  the 
2,3%  amplitude.  But  the  density  profiles  at  the  time  when  all  the 
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Figure  5.10  TQIPERATURE  PROFILE  OF  A "STABLE"  LOOP 
AT  VARIOUS  TIMES  IN  ITS  EVOLUTION. 
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plasma  in  each  loop  reaches  3 X 10  K are  very  different.  Comparing 
Figure  ^.9  and  5.11,  we  note  that  the  density  profile,  for  the  loop 
with  the  2,^%  perturbation  amplitude,  at  t ■ 5^0  sec  Is  much  smoother 
than  the  profile  at  the  corresponding  time,  t » 570  sec.  In  the  case 
of  the  5^  amplitude.  For  the  loop  with  the  smaller  amplitude,  the 
density  at  the  top  Is  less  than  twice  that  at  the  sides  and  the 
compressed  region  Is  of  much  larger  size  scale.  This  loop  would  not 
appear  to  have  a very  bright,  compact  knot  In  Ha;  It  would  appear  to 
be  of  approximately  uniform  brightness  Instead.  It  would  correspond 
to  a loop  prominence  event  that  does  not  begin  as  a bright  knot. 

For  Initial  amplitudes  larger  than  5^>  evolution  of  the  loop 

does  not  change  qualitatively  from  the  5^  case.  The  knot  appears 
brighter,  smaller  and  lasts  longer,  but  the  overall  behavior  Is  the 
same.  For  example,  assuming  a perturbation  amplitude  at  Instability 
onset  of  10^  In  a loop  with  the  same  parameters  as  before,  results  In 
relatively  larger  density  ratios  between  the  compressed  region  and 
the  sides,  (>  50  Instead  of  IG).  In  this  case  the  temperature  variations 
last  for  180  sec,  which  Is  a significant  fraction  of  the  lifetime  of 
knots. 

The  maximum  velocities  in  the  lOf?  case  occur  at  the  condensation 
boundary  as  in  the  previous  cases.  They  attain  a value  of  -1.4X10^ 
cm/ sec,  implying  a Mach  number  of  M = 4.  For  the  case  of  a ampli- 
tude, we  obtained  M = 2,  and  for  the  2.55^  case,  the  maximum  velocity 
is  3»2  X 10^  cm/ sec,  i.e.  M 1.  Hence,  the  velocity  is  approximately 
proportional  to  the  perturbation  amplitude  (for  small  amplitudes)  as 
expected  from  the  analytic  model,  equation  (5.48). 
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5.5  Post-Cooling  Stage 

We  define  the  third  stage  of  the  evolution  of  a oost-flsre  loop 

to  be  that  phase  during  which  the  temperature  remains  approximately 
4 

constant  at  3 X 10  K,  Section  5.1*  This  evolution  Is,  of  course, 

determined  by  the  pressure  and  velocity  profiles  that  are  present  at 

the  end  of  the  radiative  cooling  stage  and  by  gravity.  We  find  that 

there  are  two  possible  general  forms  that  the  third  stage  evolution  can 

take,  depending  on  whether  thermal  Instabilities  have  occurred. 

If  the  loop  Is  stable,  l.e.  no  condensations  form,  then  as  noted 

In  the  previous  section  the  density  profile  Is  approximately  uniform 

at  the  time  that  the  cooling  ends.  In  addition,  the  velocities  generated 

by  the  cooling  process  are  small,  « c.  Hence,  the  subsequent  behavior 

of  the  plasma  Is  dominated  by  gravity.  Since  the  gravitational  scale 
kT 

height,  , at  chromospheric  temperatures  Is  much  less  than  the  size 

%% 

of  typical  loops,  the  loop  plasma  begins  to  fall  freely.  The  exact  motion 
Is  complicated  due  to  the  variation  of  the  gravitational  force  along  the 
length  of  the  loop;  however,  as  one  would  expect,  we  find  that  the 
maximum  velocities  occur  at  the  base  of  the  loop.  They  are  of  the  order 

of  the  free-fall  velocity,  and  the  time  scale  for  their  growth  Is  of 

9 

the  order  of  the  free-fall  time.  For  a loop  of  height,  H = 3 X 10  cm, 
the  free-fall  velocity  la  v^  = = 10^  cm/ sec,  and  the  free- 
fall  time  Is  H/vj  « 300  sec. 

If,  on  the  other  hand,  the  loop  Is  unstable  and  a condensation 
forms,  then  we  find  that  large  pressure  gradients  are  created  which 
have  an  Important  effect  on  the  velocity  profile.  Supersonic  velocities 
form  near  the  top  of  the  loop  well  before  gravitational  effects  become 
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significant,  and  at  heights  where  the  free-fall  velocity  is  small.  They 
are  a result  of  the  rapid  expansion  of  the  compressed  region.  In  this 
section  we  describe  in  detail  the  post-cooling  evolution  of  the  unstable 


loop  discussed  in  Section  5«^.3  and  compare  it  with  the  evolution  cf 
a stable  loop. 

The  radiative  cooling  stage  ends  at  t = 570  sec  for  the  unstable 
loop  with  the  5'^  perturbation  amplitude.  Section  5*^«3.  '*'he  density 

profile  at  this  time  is  shown  in  Figure  5»9  velocity  profile 

in  Figure  5*12.  The  velocities  are  residual  from  the  non-uniform 
cooling  stage;  hence,  within  25,000  km  of  the  top  they  are  directed 
up  the  loop.  Near  the  base  they  are  directed  downward  due  to  material 
having  condensed  onto  the  chromosphere.  The  velocity  profile  indicates 
that  the  condensation  continues  to  be  compressed  (and  brighten  in  Ha) 
for  some  time  after  the  cooling  ends. 

Figure  5*9  implies  the  existence  of  large  pressure  gradients  which 
act  to  reverse  the  upward  velocities.  Since  the  temperature  is  uniform, 
the  pressure  profile  is,  of  course,  directly  proportional  to  the  density. 
The  time  required  to  reverse  the  velocity  is  o£  order 


pyX 

" P * 


(5.53) 


.8 


where  X is  the  width  of  the  compressed  region,  S lo  cm  in  our  case. 

,6 


From  Figure  5.12,  the  upward  velocity  is  approximately  - 6 X 10  cm/sec; 
hence,  M « 2 and  equation  (5.53)  yields,  T ■ 60  sec.  This  time  is 
shorter  than  the  time  scale  over  which  gravity  is  effective. 

In  agreement  with  this  estimate  for  t,  we  find  that  the  velocities 
reverse  direction  at  t - 620  sec  (50  sec  after  the  cooling  ends).  The 


87 


V (cm /sec) 


30  40 

LOOP  BASE 


LOOP  TOP 


Figure  5.12  VELOCITY  PROFILES  OF  A "STABLE" 

AND  AN  "UNSTABLE"  LOOP  AT  VARIOUS  TIMES  IN 
THEIR  EVOLUTION.  The  broken  lines  refer 
to  the  velocity  profile  of  the  unstable  loop 
at  the  onset  of  the  post-cooling  stage,  t * 57' 
sec.  The  broken  line  with  no  dots  indicates 
negative  velocities,  directed  ^ the  loop,  the 
line  with  dots  indicates  positive  velocities, 
directed  down  the  loop. 
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compression  reaches  its  maximum  density,  6 X 10  cm  , and  its  minimum 
width,  < ^ X 10^  cm,  at  this  time.  Figure  ^.13«  We  note  that  the 
minimum  density  in  the  loop  is  approximately  2 X 10^^  cm  ^ at  a position 

Q 

1.3  X 10  cm  from  the  top.  Therefore,  the  brightness  in  Hat  of  the 
compressed  region  should  be  almost  3 orders  of  magnitude  larger  than 
that  of  the  sides  of  the  loop;  however,  since  the  condensation  is 
optically  thick  in  HOf,  radiative  transfer  effects  limit  its  brightness. 

We  also  plot  the  density  profile  at  t = 620  sec  as  a function  of 
Lagrangean  coordinate,  x.  Figure  ^.l4.  Note  that  this  profile  is 
much  smoother  than  the  corresponding  one  in  Figure  5.13  ~ the  compressed 
region  is  clearly  visible  and  is  of  almost  uniform  density.  Over 
30^  of  the  loop  plasma  is  in  the  condensation;  however,  it  occupies 
leas  than  10  of  the  loop  volume.  Figures  5*13  and  5«l4  clearly 
Illustrate  the  necessity  of  using  Lagrangean  coordinates  for  our 
computer  calculations.  If  we  were  to  attempt  to  follow  the  density 

7 

in  Eulerian  coordinates,  we  would  need  a grid  spacing  finer  than  10  cm, 

■i 

which  implies  using  almost  10  points  over  the  length  of  the  loop.  In 
fact,  we  only  use  64  points  with  the  Lagrangean  grid. 

The  density  profiles  of  the  loop  at  t ■ 620,  720,  840  and  960  sec 
are  plotted  in  Figure  5«13«  The  rapid  expansion  of  the  condensation 
can  clearly  be  seen;  the  velocity  of  the  condensation  boundary  is 
approximately  5 X 10^  cm/sec  which  is  supersonic. 

We  plot  the  velocity  profile  at  t ■ 840  in  Figure  5>12  and,  for 
comparison,  we  also  plot  the  profile  at  t > 34o  sec  for  a loop  with 
the  same  geometry  but  with  no  initial  perturbation,  l.e.  a completely 
stable  loop.  Below  10^  cm  from  the  top,  the  velocities  in  the  stable 
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Figure  ^.13  density  PROFILE  OF  AN  "UNSTABLE"  LOOP 
AT  VARIOUS  TIMES  IN  THE  POST-COOLING  PHASE  OF 


ITS  EVOLUTION. 
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Figure  5,14  DENSITY  IN  AN  "UNSTABLE"  LOOP  AS 
A FUNCTION  OF  LAGRANGEAN  COORDINATE,  x, 

AT  THE  TIME  WHEN  THE  CONDENSATION  DENSITY 
IS  MAXIMUM. 


loop  are  larger  than  the  ones  in  the  unstable  case.  This  Is  because 
no  upward  velocities  are  generated  during  the  radiative  cooling  stage 

Q 

In  the  stable  case.  On  the  other  hand,  within  10  cm  of  the  top,  the 
velocities  In  the  unstable  loop  are  higher  due  to  the  expansion.  The 

g 

maximum  expansion  velocicy  at  s « 8 X 10^  cm  Is  almost  equal  to  the 
maximum  free-fall  velocity  at  s = 32  X lo”  cm  for  the  unstable  loop. 

If  viewed  In  an  Htt  movie,  these  two  loops  would  appear  to  have  quite 
a different  behavior.  In  the  stable  loop  material  would  appear  to 
accelerate  uniformly  as  It  fell  down  the  sides,  while  In  the  unstable 
loop,  the  plasma  would  seem  to  explode  out  of  the  HOT  bright  knot, 
reaching  a large  fraction  of  Its  maximum  velocity  near  the  top. 

Once  the  condensation  expands  sufficiently  that  pressure  gradients 
become  negligible,  the  subsequent  evolution  of  the  plasma  is  dominated 
by  gravity.  For  the  unstable  loop,  we  find  that  the  pressure  (and 
hence  the  density)  becomes  iMnotonlcally  decreasing  with  height  by 
t * 20  min.  The  density  at  the  top  Is  only  1.9  X 10^^  cm  ^ at  this 
time.  By  t ■ 23  min,  the  top  density  falls  below  10^^  cm  ^ and  we 
terminate  our  calculation.  For  the  stable  loop,  this  density  Is  reached 
sooner  at  t - 18  min;  however,  the  maximum  free-fall  velocity  in  both 
cases  Is  approximately  10  cm/sec. 

In  summary,  we  find  that  the  main  difference  In  the  evolution  of 
a stable  and  unstable  loop  Is  In  the  appearance  of  bright  knots  and  In 

the  acceleration  of  cold  plasma  near  the  top.  Since  we  somewhat 

4 

arbitrarily  fixed  the  temperature  at  3 X 10  K,  the  details  of  our 
results  on  the  post-cooling  phase  of  flare  loops  may  be  Inaccurate, 
but  we  believe  that  the  Important  features  are  correct.  We  compare 
these  features  with  observations  In  the  following  section. 
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6.  DISCUSSION 
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6. 1 Comparison  with  Observations 

If  we  compare  the  results  In  the  previous  chapter  on  the  evolution 
of  an  unstable  loop  with  the  description  of  loop  prominences  given  In 
chapter  2,  we  note  that  they  are  In  good  agreement.  For  the  case  of  the 
loop  with  a Initial  perturbation  amplitude,  the  average  density  In  the 

11  _3 

condensation  during  Its  lifetime  Is  ^ 10  cm  , and  Its  size  Is  small 

0 

compared  to  the  size  of  the  loop,  < 5 X 10°  cm.  In  Hof  It  would  appear 
as  a bright  knot.  The  plasma  velocities  near  the  condensation  are  large; 
there  Is  a shock  at  the  boundary  during  Its  formation  and  a large 
expansion  velocity  during  Its  decay.  This  agrees  with  Jefferies  and 
Orrall's  (196^a)  observation  of  large  streaming  velocities  at  the  top 
of  active  loop  prominences.  Defining  the  lifetime  of  a knot  as  the 
time  Interval  from  the  Instant  It  first  appears  in  Hv  to  when  It  expands 
to  a size  of  half  the  loop  length;  we  find  that  the  lifetime  In  our 
case  Is  approximately  8 minutes.  Bruzek  and  Kuperus  (1972)  quote  the 
lifetime  of  a knot  as  3 - 10  minutes. 

The  maximum  velocity  that  we  obtain  at  the  loop  base,  113  km/ sec, 
equals  the  free-fall  velocity.  The  average  velocity  Is  slightly  lower, 

< 100  km/sec.  In  agreement  with  observations.  There  Is  a shock  at  the 
base  where  the  falling  matter  Impacts  with  the  chromosphere;  however, 
the  radiative  losses  are  sufficient  to  dissipate  the  mechanical  energy 
flux.  In  an  unstable  loop,  large  velocities  appear  near  the  top  (at 
the  knot  boundary)  due  to  pressure  gradients;  whereas.  In  a stable  loop 
material  accelerates  uniformly  down  the  sides.  We  calculate  the  total 
lifetime  of  a loop  event  as  approximately  30  minutes.  Bruzek  (196^a) 
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gives  the  observed  lifetime  of  typical  loops  as  30  - 60  minutes.  Of 
course,  these  time  scales  and  velocities  viU  vary  depending  on  the 
size  and  density  of  a loop,  equations  (5*2)  and  (5.3).  lower 

densities  and  larger  loop  heights,  the  radiative  cooling  time  and 
the  free-fall  time  increase.  However,  the  long  lifetime  of  loop 
prominence  systems,  ^ 10^  sec,  does  not  appear  to  be  due  to  a cooling 
effect.  Even  if  all  conductive  losses  are  suppressed  by  a large 
magnetic  field  compression  factor  T,  the  radiative  losses  alone  for 
densities  of  ~ 10^^  cm"^  and  temperatures  of  ~ 10^  K imply  a cooling 
time  of  only  10^  sec.  Hence,  the  long  lifetime  of  a loop  system  is 
most  likely  due  to  a long-lived  flare  heating  process. 

In  general  we  find  that  the  temperature,  densities  and  velocities 
measured  in  a loop  prominence  are  accurately  reproduced  by  our  model. 

From  this,  we  conclude  that  the  thermal  instability  mechanism  by  itself 
is  sufficient  to  account  for  a loop  prominence  event.  The  appearance  of 
knots  and  their  subsequent  evolution  can  be  understood  as  being  due  to 
the  temperature  and  density  dependence  of  the  cooling  function  (Figure  1.1) 
for  the  coronal  plasma.  Our  only  important  assumption  is  the  form  of 
the  temperature  and  density  profiles  during  the  radiative  cooling  stage. 

If  a perturbation  of  small  but  finite  amplitude  is  present  a condensation 
will  form.  If,  on  the  other  hand,  the  profiles  are  uniform  as  we  would 
expect  from  the  discussion  in  Section  5.2.1,  then  no  knots  appear  and 
the  evolution  of  the  loop  plasma  will  more  closely  resemble  that  of 

coronal  rain. 

We  have  given  no  Justification  for  choosing  the  particular  form 
of  the  perturbation  in  equation  (5.‘+9),  e««Pt  that  it  is  the  one  with 


the  largest  wavelength  and,  hence,  the  fastest  growth  rate.  In  accordance 
with  observations,  condensations  may  form  at  the  sides  of  loops  If  the 
perturbation  profiles  are  suitably  chosen;  however,  these  knots  are  not 
as  dense  or  as  long-lived  as  those  centered  at  the  top.  Unfortunately, 
the  observations  are  not  accurate  enough  to  determine  the  temperature 
and  density  profiles  of  a loop  to  within  hence  we  cannot  directly 
check  our  assumption.  The  soft  x-ray  emission  Is  observed  to  be 
greater  at  the  tops  of  loops  (Cheng  and  Wldlng,  1975;  Pallavaclnl  et  al., 
1975 )»  but  that  would  be  expected  due  to  the  area  factor  even  If  the 
temperature  were  slightly  lower  there. 

There  Is  reason  to  expect  that  the  temperature  will  be  higher  at 
the  sides  of  a flare  loop  If  the  heating  Is  due  to  high  energy  electrons 
as  In  the  models  proposed  by  Sturrock  (I968),  Petrosian  (1973)  and  others 
(Svestka,  1976,  Section  3>2).  In  these  models  the  bulk  of  the  flare 
energy  Is  carried  down  the  field  lines  by  non-thermal  electrons  and 
deposited  In  the  chromosphere  or  lower  corona.  Hence,  the  temperature 
should  Initially  be  higher  near  the  base  of  a loop.  If  the  electrons 
have  non-zero  pitch  angles  and  are  trapped  In  the  flux  tube,  then 
again,  we  would  expect  higher  temperatures  at  the  sides  of  a loop.  The 
electrons  spend  more  time  near  their  mirror  points  than  near  the  top  of 
a loop,  since  their  velocity  parallel  to  the  field  Is  highest  at  the  top. 
This  process  also  has  the  attraction  that  If  a loop  Is  symmetrical  about 
the  top,  then  the  heating  due  to  trapped  energetic  electrons  will  have 
the  same  symmetry;  hence  any  temperature  perturbation  will  be  symmet- 
rical about  the  top  In  agreement  with  our  assumption  In  equation  (5«^9)> 
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intend  to  modify  the  boundary  conditions  to  permit  an  upward  mass  flux 
at  the  base  of  the  loop.  This  will  allow  us  to  study  the  evaporation 
of  chromospheric  material,  and,  perhaps,  to  determine  the  form  of 
flare  heating. 

6.2  Conclusions 

To  summarize,  we  restate  the  main  conclusions  of  our  calculations 

here: 

1.  For  typical  post-flare  loops,  the  conductive  damping  time  at  high 

7 

temperature,  ^10  K,  is  so  short  that  all  temperature  and  pressure 
gradients  are  damped  quickly  compared  to  the  cooling  time.  Hence  the 
analytic  models  proposed  by  Antiochos  and  Sturrock  (19?6a,  b)  adequately 
represent  the  initial  cooling. 

o ' 

2.  The  Inclusion  of  the  area  factor  does  not  alter  the  linear  growth 
rate  of  thermal  perturbations,  equation  (l,l4). 

3.  The  effect  of  the  nonlinearities  is  to  increase  the  stability  of 
the  plasma.  Figure  5.3. 

U.  Observable  condensations  will  grow  if  criteria  1 - 4 in  Section  ^.3 
are  satisfied.  The  important  requirement  is  that  the  perturbation 
amplitude  at  the  onset  of  instability  be  ^ 

3m  Shocks  form  at  the  boundary  of  knots  during  their  compression  phase. 

Large  velocities  are  generated  by  their  subsequent  expansion. 

6.  The  appearance  of  bright  knots  is  mostly  due  to  density  variations  ] 

along  the  loop,  since  temperature  variations  last  for  only  a short  time, 

< 100  sec.  *: 

7.  The  evolution  of  loop  prominences  can  be  understood  as  simply  due  ^ 

to  thermal  instability  in  chromospheric  material  that  has  been 
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across  field  lines  are  not  required. 

8.  The  radiative  cooling  time  for  temperatures  and  densities  typical  of 
post-flare  loops  is  too  short  to  account  for  the  long  lifetime  of  loop 
prominence  systems;  hence,  a long-lived  (time  scales  of  the  order  of 
several  hours),  flare  heating  process  is  required. 


APPENDIX  A 
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FLUID  EQUATIONS  FOR  A DIPOLE  LOOP 


We  wish  to  derive  the  fluid  equations  for  the  case  of  a flux  tube. 
Since  we  assume  that  the  field  dominates  the  plasma,  we  can  most  easily 
Incorporate  the  effects  of  the  field  by  letting  the  field  lines  define 
a curvilinear  coordinate  system.  This  Is  always  possible  If  the  field 
Is  potential  (Morse  and  Feschback,  1953t  p. l4). 


12  3 1 

Let  (x  , X , x^)  denote  the  coordinate  system,  where  x » s,  the 


distance  along  the  field.  We  are  only  Interested  In  an  Individual 
loop  whose  cross-sectional  area  Is  very  small  compared  to  the  radius 
of  curvature  of  the  field.  Thus,  we  can  approximate  the  metric  In 
the  loop  by; 


0 

0 


(Al) 


A(x^) 


To  obtain  the  fluid  equations  In  terms  of  our  coordinate  system,  we  first 
rewrite  them  In  covariant  form.  They  are; 


fj+(nv^)^  = 0 . 


(A2) 


+ (nmjjv^v^  + p6|  + . nm^g^  , 

ail  i p)  + I pv^  + qJ  + 


(A3) 


at 


ntUjjg^v  + J- 


where  we  have  used  the  form  of  the  equations  given  by  Braglnskll  (I965) 
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H 


(A4); 


and  the  tensor  notation  of  Adler,  Bazin  and  Schlffer  (I965),  Chapter  2. 
Note  that  all  spatial  derivatives  are  covariant  derivatives  and  that  the 
Einstein  convention  of  summing  over  repeated  Indices  Is  used. 


I 


I 


The  equations  A2  - A4  simplify  considerably  In  our  case  since  we  have 


2 3 

no  variables  dependent  on  x or  x and  all  vectors  have  a component  in 


the  x^  direction  only.  Hence, 

v^  = (v  (x^),  0,  0)  . 


(A5) 


The  heat  flux  and  the  stress  tensor  are  given  by  Braginskii 


(1965); 


■ (■'  °) 


(A6) 


and 


""ij  = tniJ  ^ "Jlli  ■ 3 *ij  "„k) 


(A?) 


hence. 


^ _ 2 1 ^(Av) 


dx 


.(A8) 


Using  equation  (A^)  and  the  metric  (Al),  the  mass  equation  (A2)  becomes: 

(A9) 


^ 4.  i ^Anv)  ^ 

dt  A as  “ ^ * 


where  we  have  replaced  x^  by  s.  Using  (A9)  and  (a8),  the  momentum 


equation  (A3)  becomes: 


99 


> I 


(AlO) 


w 


2dv  _ V 3 £a  A 
3s  3s 


(All) 


and  we  have  replaced  g,  by  g . 

1 II 

Finally,  using  equations  (AlO),  (A9)  and  (a6),  the  heat  equation 
(aU)  can  be  expressed  as: 

2 


^(Av)  _ 1 ^(1^)  Jlw'^ 
2 \ 3t  3s  / 2 A 3s  A 3s  3s  3 


f . (A12) 


Since  T|  is  small,  we  can  neglect  the  viscosity  terms  to  obtain  the  set 
of  equations  (3.20)  - (3.22)  in  Section  3,6. 

We  now  need  to  calculate  A(s)  and  8(|(®^  ® dipole  field.  Using 

the  coordinates  defined  in  Figure  (3.1),  a dipole  field  is  given  by 
(Jackson,  I962,  p,  143), 


B a (2  sin  6 2^  - cos  e tg) 


(AI3) 


The  equation  of  a field  line  of  height  R is  given  by: 

p 

r/R  = cos  0 . 

Letting 

V = sin  6 , 

the  distance  along  the  field  line, 

di  + r^(de)^  , 
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(A14) 


(A15) 


(Al6) 


t 


Figure  A.l  PROFILES  OF  THE  AREA  OF  A LOOP  DUE  TO  A 
DIPOLE  SOURCE  AND  OF  THE  COMPONENT  OF  GRAVITY  PARALLEL 
TO  THE  LOOP.  The  solid  line  refers  to  A(S)  end  the 
I _ broken  line  to  g.<(S).  Both  A end  gj.  ere  noraellsed  to 

^ ‘ equel  unity  at  tnelr  respective  mexime. 
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APPENDIX  B 
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C B.l  FINITE  DIFFERENCE  EQUATIONS 
C 

C DETAILS  OF  DIFFERENCING  SCHEME  FOR  EQUATIONS  (4.3  - 4.7)  WITH 
C BOUNDARY  CONDITIONS  (4.18-4.19).  THE  FOLLOWING  IS  THE  CODE  TO 
C ADVANCE  VALUES  OF  VARIABLES  BY  ONE  TIME  STEP,  DT.  THE  VALUE  OP  DT 
C FOR  THE  NEXT  TIME  STEP  CAN  BE  EVALUATED  USING  CONDITION  (4.4o). 

C THAT  ROUTINE  IS  NOT  INCLUDED  HERE.  IN  THE  NEXT  SECTION  A ROUTINE 
C FOR  TABLE  LOOK-UPS  IS  GIVEN. 

C 

C WE  DIVIDE  EACH  LEG  OF  THE  LOOP  INTO  32  EQUAL  MASS  SLABS.  EACH  SLAB 
C IS  OF  WIDTH  (MASS)  DX.  WE  CHOOSE  THE  POINT  X - 0 TO  BE  THE  TOP  OF 
C THE  LOOP. 

C 

C LIST  OF  VARIABLES  USED  IN  CODE: 

C ALL  ARRAYS  HAVE  33  ELEMENTS.  (SOMETIMES  NOT  ALL  ELEMENTS  USED). 

C PRIMARY  VARIABLES : 

C VALUES  OF  PLASMA  PARAMETERS  AT  TIME  » T.  NEEDED  AS  INPUT  TO  ADVANCE 
C TIME  STEP.  MUST  BE  OUTPUTTED  TO  USE  FOR  THE  NEXT  STEP. 

C "S(J)"  = EULERIAN  POSITION  OF  LAGRANGEAN  POINT  (j-l)*DX 
C "V(J)"  = VELOCITY  AT  LAGRANGEAN  POINT  (j-l)*DX. 

C "T(J)"  * TEMPERATURE  AT  LAGRANGEAN  POINT  (J-1/2)*DX. 

C (NOTE  THAT  WE  POSITION  SOME  VARIABLES  DIFFERENTLY  FOR  PROPER 
C CENTERING  OF  DIFFERENCE  EQUATIONS.) 

C "P(J)"  = PRESSURE  AT  LAGRANGEAN  POINT  (j-l/2)*DX. 

C ”AN(J)"  = NUMBER  DENSITY  AT  LAGRANGEAN  POINT  (J-1/2)*DX. 

C "Q(J)"  - ARTIFICIAL  VISCOSITY  AT  LAGRANGEAN  POINT  (J-1/2)*DX. 

C 

C PROVISIONAL  VARIABLES: 

C VALUES  OF  PLASMA  PARAMETERS  AT  TIME  = T + DT/2.  THE  SPATIAL 
C POSITIONING  IS  SAME  AS  ABLVE. 

C "S1(J)",  "V1(J)",  "T1(J)",  "P1(J)",  "ANI(J)",  "Q1(J)" 

C 

C SECONDARY  VARIABLES : 

C VALUES  OF  SECONDARY  PARAMETERS  AT  EITHER  TIME  - T OR  T + DT/2. 

C MAY  BE  OUTPUTTED  IF  DESIRED. 

C "VDER(J)"  = D(V(J))/DT. 

C "TDER(J)"  . D(T(J))/DT. 

C "AS(JV'  - AREA  FACTOR  AT  LAGRANGEAN  POINT  (j-l)*DX. 

C "GS(J)"  * COMPONENT  OF  GRAVITY  PARALLEL  TO 
C LOOP  AT  LAGRANGEAN  POINT  (J-1)*DX. 

C "T12(J)"  - SQUARE  ROOT  OF  T(J). 

C "T72(J)"  - T(J)  RAISED  TO  THE  POWER  OF  7/2. 

C "ALT(J)"  - COX  AND  TUCKER  RADIATIVE  LOSSES  AT  TEMPERATURE  T(J). 

C "HFLX(J)"  - HEAT  FLUX  AT  LAGRANGEAN  POINT  (j-l)*DX. 

C 

C ARRAYS  USED  FOR  STORAGE  ONLY:  "V12(j)",  "AV(j)”,  ”D(j)". 

C 

C CONSTANTS:  CHOSEN  TO  MAKE  VARIABLES  DIMENSIONLESS. 

C "DX"  - 2/3,  "BQVl",  "QOPO",  "BQT3",  "EQT32". 

C "AIB"  - VOLUME  OF  LEG  OF  LOOP. 

C "DT"  - TIME  INCREMENT,  CONSTANT  ONLY  FOR  EACH  TIME  ADVANCE.  I 

C "DT2"  - DT/2,  "DT3"  » DT/3,  "DT4"  - DT/4.  j 
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C VARIABLES  USED  FOR  STORAGE  ONLY:  "AIL,  AIF,  DELTV". 

C 

C BOUNDARY  CONDITIONS:  VALUES  CONSTANT  FOR  ALL  TIME 
C ROUTINE  TO  INCREMENT  TIME  BY  DT/2. 

C CALCULATE  PROVISIONAL  VARIABLES:  Vl(2),  Sl(2),  ANl(l),  Ql(l) 

VDER(2)  - BQV1»AS(2)*(P(2)  + Q(2)  - P(l)  - Q(l))  + GS(2) 

VI (2)  - V(2)  + DT2*VDER(2) 

V12(2)  » Vl(2)  + V(2) 

Sl(2)  - S(2)  + DT4*V12(2) 

C— INSERT  ROUTINE  TO  CALCULATE  AS(2),  GS(2),  AIL  - SEE  SECTION  B.2 
C AIL  - VOLUME  OF  LOOP  FROM  TOP  TO  POINT  SI (2) 

310  ANl(l)  . DX/(A1L) 

AV(2)  . AS(2)*V1(2) 

DELTV  * Vl(2)  - Vl(l) 

IF (DELTV )330,  320,  320 
320  Ql(l)  « 0. 

GO  TO  340 

330  Ql(l)  - Q0P0^^AN1(1)*DELTV»DELTV 
34o  continue 

C START  MAIN  DO  LOOP  lllllllllllllllllllllllllllllllllllllll 

DO  46o  J » 1,  30 

C CALCULATE  Vl(j+2),  Sl(j+2),  ANl(jfl),  Ql(J+l) 

VDER(J+2)  - EQVl*AS(J+2)*(P(J+2)  + Q(J+2)  - P(J+1)  - Q(J+1)) 

1 + GS(j+2) 

Vl(J+2)  » V(J+2)  + DT2*VDER(j+2) 

VlE(J+2)  = Vl(J+2)  + V(J+2) 

Sl(J+2)  - S(J+2)  + DT4*V12(J+2) 

C-— INSERT  ROUTINE  TO  CALCULATE  AS(J+2),  GS(j+2),  AIF  - SEE  SECTION  B.2 
C AIF  = VOLUME  OF  LOOP  FROM  TOP  TO  POINT  Sl(j+2) 

370  AN1(J+1)  - DX/(AIF-AIL) 

AIL  * AIF 

AV(J+2)  - AS(j+2)*Vl(J+2) 

DELTV  * Vl(J+2)  - V1(J+1) 

IF (DELTV)  390,  380,  380 
380  Q1(J+1)  = 0. 

GO  TO  400 

390  A1(J+1)  = Q0P0*AN1(J+1)*DELTV*DELTV 

400  CONTINUE 
C CALCULATE  T1(J),  Pl(J) 

D(J+1)  » AS(J+1)»AS(J+1)*AN1(J)*AN1(J+1)/(AN1(J)  + ANl(j+l)) 
HFLX(J+1)  - D(J+1)*(T72(J+1)  - T72(J)) 

TDER(J)  . (P(J)  + Q1(J))*(AV(J)  - AV(J+1))  + ANl(j)*ALT(J) 

1 + BQT3*(HFLX(J+1)  - HFLX(J)) 

T1(J)  » T(J)  + DT2*TDER(J) 

IF(T1(J)  .LT.  T(33))GO  TO  44o 

C— INSERT  ROUTINE  TO  CALCULATE  T12(j),  T72(J),  ALT(j)  - SEE  SECTION  B.2 
GO  TO  450 

440  T1(J)  « T(33) 

ALT(J)  - ALT03) 

T72(J)  - T72(33) 

450  P1(J)  - AN1(J)»(T1(J)) 

460  CONTINUE 

C END  OF  MAIN  DO  LOOP  111111111111111111111111111111111111111 
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C CALCULATE  ANl(32),  Q1C32) 

AN1(32)  - DX/(AIB-AIL) 

DELTV  = VI (33)  - VI (32) 

IF(DELTV)  480,  1+70,  1+70 
1+70  Qlb2)  - 0. 

GO  TO  1+90 

l+80  Ql(32)  - Q0P0*AN1  (32  )*DELTV»DELTV 
490  CONTINUE 
C CALCULATE  Tl(3l),  PI (31) 

D(32)  - AS(32)*AS(32)*AN1(31)*AN1(32)/(AN1(32)  +ANl(3l)) 

HFLX(32)  - D(32)*(T72(32)  - T72(31)) 

TDER(31)  = (P(31)  + Q1(31))*(AV(31)  - AV(32))  + ANl(31)*ALT(31) 

1 + EQT3*(HFLX(32)  - HFLX(3l)) 

Tl(31)  = T(31)  + DT2*TDER(31) 

IF(T1(31)  .LT.  T(33 ))G0  TO  530 

C— INSERT  ROUTINE  TO  CALCULATE  T12(31),  T72(31),  ALT(31)  - SEE  SECTION  B.2 
GO  TO  54o 

530  T1(31)  = T(33) 

ALT(31)  = ALT(33) 

T72(31)  = T72(33) 

540  PI (31)  - AN1(31)*(T1(31)) 

C CALCULATE  Tl(32),  PI (32) 

hflx(33)  - AS (33)^(172(33)  - tt2(32))/(s(33)  - SI (32)) 

TDER(32)  . (P(32)  + Q1(32))*AV(32)  + ANl (32)*ALT(32) 

1 + EQT32*HFLX(33)  - EQT3*HFLX(32 ) 

Tl(32)  = T(32)  + DT2*TDER(32) 

IF(T1(32)  .LT.  T(33))  go  to  58O 

C— INSERT  ROUTINE  TO  CALCULATE  T12(32),  T72(32),  ALT (32)  - SEE  SECTION  B.2 
GO  TO  590 

580  Tl(32)  = T(33) 

ALT(32)  - ALT(33) 

T72(32)  = T72(33) 

590  P1(32)  = AN1(32)»(T1(32)) 


C END  OF  TIME  INCREMENT  BY  DT/2.  BEGIN  INCREMENT  BY  DT. 

C CALCULATE  PROVISIONAL  VARIABLES:  V(2),  S(2),  AN(1),  Q(l) 

VDER(2)  - EQV1*AS(2)*(P1(2)  + Ql(2)  - Pl(l)  - Ql(l))  + GS(2) 

V(2)  - V(2)  + DT*VDER(2) 

S(2)  - S(2)  + (V12(2)  + V(2))*DT3 

C— INSERT  ROUTINE  TO  CALCULATE  AS(2),  GS(2),AIL  - SEE  SECTION  B.2 
C AIL  « VOLUME  OF  LOOP  FROM  TOP  TO  POINT  S(2) 

620  AN(1)  - DX/(AIL) 

DELTV  - V(2)  - V(l) 

IF (DELTV )640,  630,630 
630  Q(l)  = 0. 

GO  TO  650 

640  Q(l)  - Q0PC*AN(1)^^DELTV*DELTV 
650  CONTINUE 

C START  OF  MAIN  DO  LOOP  2222222222222222222222222222222222222 
DO  770  J-1,  30 

VDER(J+2)  « EQVl*AS(J+2)*(Pl(J+2)  + Ql(J+2)  - Pl(J+l)  - Ql(J+l)) 

1 + GS(J+2) 

V(J+2)  - V(J+2)  + DT*VDER(j+2) 

S(J+2)  - S(j+2)  + (V12(J+2)  + V(J+2))*DT3 
C— INSERT  ROUTINE  TO  CALCULATE  AS(J+2),  GS(J+2),  AIF  - SEE  SECTION  B.2 
C AIF  - VOLUME  OF  LOOP  FROM  TOP  TO  POINT  S(J+2) 


680  AN(J+1)  » DX/(AIF-AIL) 

AIL  - AIF 

DELTV  = V(j+2)  - V(J+1) 

IF(DELTV)700,  690,  690 
690  Q(J+1)  * 0. 

GO  TO  710 

700  Q(J+1)  « QOPO^AN  ( J+1  )*DEL’mDELTV 
710  CONTINUE 
C CALCULATE  T(j),  P(J) 

HFLX(J+1)  - D(J+1)*(T72(J+1)  - T72(J)) 

TDER(J)  - (P1(J)  + Q1(J))*(AV(J)  - AV(J+l))  + ANl (j)*ALT(j) 

1 + EQT3*(HFLX(J+1)  - HFLX(J)) 

T(J)  - T(J)  + DT*TDER(J) 

IF(T(J)  .LT,  T(33))GO  TO  750 

C— INSERT  ROUTINE  TO  CALCULATE  T12(J),  T72(J),  ALT(j)  - SEE  SECTION  B.2 
GO  TO  760 
750  T(J)  - T(33) 

T12(J)  - T12(33) 

T72(j)  = T72(33) 

760  P(J)  - AN(J)*(T(J)) 

770  CONTINUE 

C END  OF  MAIN  DO  LOOP  2222222222222222222222222222222222222222 
C CALCULATE  AN(32),  Q(32) 

AN (32)  = DX/(AIB-AIL) 

C CALCULATE  T(31),  P(3l) 

DELTV  = V(33)  - y{?,2) 

IF (DELTV ) 790 ,780,780 
780  Q(32)  = 0 

GO  TO  800 

790  Q(32)  » Q0P0»AN(32)*DELTV*DELTV 
800  CONTINUE 

HFLX(32)  - D(32)»(T72(32)  - T72(31)) 

TDER(31)  - (P1(31)  + Ql(3l))*(AV(31)  - AV(32))  + AN1(31)*ALT(31) 

I + EQT3*(HFLX(32;  - HFLX(3l)) 

T(31)  - T(31)  + DT»TDER(31) 

IF(T(31)  .LT.  T(33))GO  TO  840 

C-— INSERT  ROUTINE  TO  CALCULATE  T12(31),  T72(31),  ALT(3l)  - SEE  SECTION  B.2 
GO  TO  850 
8U0  T(31)  - T(33) 

ALT(31)  - ALT(33) 

T12(31)  = T12(33) 

T72(31)  - T72(33) 

850  P(31)  - AN(31)*(t(31)) 

C CALCULATE  T(32),  P(32) 

855  HFLX(33)  - AS(33)*(T72(33)  - T72(32))/(S(33)  - Sl(32)) 

TDER(32)  - (PI (32)  + Q1(32))*AV(32)  + AN1(32)*ALT(32) 

1 + BQT32*HFLX(33)  - BQT3*HFLX(32) 

T(32)  - T(32)  + DT»TDER(32) 

IF(T(32)  .LT.  T(33))GO  TO  89O 

C— INSERT  ROUTINE  TO  CALCULATE  T12(32),  T72(32),  ALT(32)  - SEE  SECTION  B.2 
GO  TO  900 
890  T(32)  - T(33) 

ALT(32)  - ALT(33) 

T12(32)  - T12(33) 

T72(32)  - T72(33) 

900  P(32)  - AN(32)»(T(32)) 

C END  OF  ROUTINE 

C ALL  PRIMARY  ^VARIABLES  ^^VE^  BEEN^  INCIffiMEmD  ^ Ifif  n o n n h k k kk  m k k w k k k x x < k kw  it 
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C B.2  ROUTINE  FOR  TABLE  LOOK-UPS 
C 

C AS  DISCUSSED  IN  SECTION  4.4  WE  MUST  EVALUATE  3 FUNCTIONS  OF  EULERIAN 
C POSITION  S(J);  AS(J),  GS(J),  AND  AIF  (OR  AIL),  AND  2 OF  TEMPERATURE 
C T(J):  T12(J)  AND  ALT(J).  WE  FIRST  TABULATE  THESE  5 FUNCTIONS  AT  240 
C LINEARLY  SPACED  VALUES  IN  EACH  OF  THE  FOUR  INTERVALS: 

C S(J)  OR  T(J)  . 1/4096  - 1/256,  1/256  - 1/16,1/16  - 1,  AND  1 - I6, 

C AND  THEN  USE  LINEAR  INTERPOLATION  TO  FIND  THE  VALUE  OF  A FUNCTION  FOR 
C A PARTICULAR  S(j)  OR  T(j).  BY  SUITABLY  NORMALIZING  S(J)  AND  T(j)  WE 
C CAN  USE  THE  FIRST  2 BYTES  OF  S(j)  OR  T(j)  TO  LOCATE  THE  NEAREST  ENTRY  IN 
C THE  TABLE.  FOR  EXAMPLE,  WE  LIST  BELOW  THE  ROUTINE  FOR  FINDING  ALT(j) 

C AND  T12(J).  IN  THIS  ROUTINE  THE  TABLES  ARE  IN  ARRAYS  AL  AND  TSQ,  THE 
C POSITION  OF  THE  NEAREST  ENTRY  IS  EITHER  "NST”  OR  "NST+1". 

C 

real*4  tch,  stref,  tabsep 
LOGICAL*!  QQ,QL(4) 

INTEGER*4  NFSET 
INTEGER*2  NN,MM 

EQUIVALENCE  (TCH,  NN,  QQ),  (TABSEP,  QL(l)),  (MM, STREF) 

DATA  NFSET/ZOOOO3EOF/,  TABSEP/ ZOOOIOOOO/, 

1 S TREF / ZOOOOOOOO / 

C-  — ROUTINE  TO  CALCULATE  T12(j),  T72(j),  ALT(j) 

C FIRST  WE  FIND  "NST".  THE  POSITION  OF  THE  ENTRY  NEAREST 
C TO  T(J)  AND  LESS  THAN  T(J). 

TCH  - T(J) 

NST  = NN  - NFSET 

C NCW  WE  FIND  THE  SEPARATION  BETWEEN  T(J)  AND  THE  POSITION  "NST". 

QL(1)  = QQ 
MM  - NN 

DIFST  - TCH-STREF 

C NOW  DECIDE  WHETHER  "NST"  OR  "NST+1"  IS  CLOSER  TO  T(j). 

C FOR  EXTRA  ACCURACY  WE  USE  SECOND  ORDER  INTERPOLATION  FOR  T12(J). 

IF (DIFST4DIFST-TABSEP)720,720,730 
720  ALT(J)  - AL(NST)  + DIFST* (AL (NST+1 ) - AL (NST )) /TABSEP 
TRATO  - DIFST/ (STREF+STREF) 

T12(J)  - TSQ(NST)*(1.+TRAT0*(1.-.5*TRAT0)) 

GO  TO  T40 

730  DIFST  - DIFST-TABSEP 

ALT(J)  - AL (NST+1)  + DIFST* (AL (NST+1 ) - AL (NST )) /TABSEP 
TRATO  - DIFST/ (STREF+STREF+TABSEP+TABSEP) 

T12(J)  - TSQ (NST+1 )*(1.+TRAT0*(1.-.5*TRAT0)) 

740  T72(J)  - TCH*TCH»TCH*(T12(J) 

C--END  OF  ROUTINE 
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w flare  event  have  been  pi^posed  as  a mechanism  for  these  types  of  active  ^eminences 
> rtfortMa  have  Investigate^  the  cooling  of  post-flare  plasmas,  Eind  attempts  to 
explain  the  formation  of  loop  prominences  as  due  to  a thermal  Instability,  ^t 
K temperatures  hotwaen  10*^  to-  10' -4^  the  solar  plasma  Is  known  to  be  unstable  to 
thermal  perturbations  because  the  radiation  losses  of  this  gas  Increase  with 
decreasing  temperature^  A simple  perturbation  analysis  Is  insufficient  to  account 
/for  the 'evSlutlon  ot  the  solar  plasma,  because  observed  differences  of  temperature 
and  density  exceed  an  order  of  magnitude.  Additionally,  the  plasma  Is  always  In 
a non-equlllbrlum  state,  so  that,  we  must  solve  numerically  the  equations  of 
motion  and  heat  transfy.  'is  ci<Ve(o/5f(^ 

— Wa^hav  develope^^one-dlmenslonal  model^for  active  loop  prominences. 

Magnetic  fields  present^n  post-flare  regions  are  strong  enough  to  dominate  the 
plasma;  hence,  only  motion  and  heat  fluxes  parallel  to  the  field  need  to  be 
considered.  The  relevant  size  scales  and  time  scales  are  such|^|t  single- 
fluid  MHD  equations  are  valid.  We  have  included-  in  ^e  modeiA^e  effects  of 
gravity,  the  geometry  of  the  field  and  conduction  loTles  to  the  chromosphere.  A 
computer  code  for_the  solution  of  our  set  of  equations  has  been  constructed.^ 
Basically,  we  treat  the  system/^as  an  Initial  value  problem  (with  certain  boundary 
conditions  at  the  chromosphere-corona  transition  region),  and  use  a two-step 
time  differencing  scheme. 

Our  calculations  Indicate  that  the  thermal  instability  mechanism  Is,  by 
Itself,  sufficient  to  account  for  the  behavior  of  active  loop  prominences.  Under 
certain  conditions.  Initial  perturbations  In  the  temperature  and  density  profiles 
of  small  amplitude  (<  and  large  size  scale  2 X 10^  cm.)  can  grow  Into 
condensations  with  temperature  and  density  differences  of  over  an  order  of 
magnitude  and  size  scales  of  less  than  10°  cm.  In  agreement  with  observations, 
the  conditions  that  must  be  satisfied  are  such  that  Loop  Prominence  Systems  are 
likely  to  occur  only  In  large  flares.  The  velocities,  densities,  and  lifetimes 
-4hat  we  obtain  for  the  loop  material  are  also  In  agreement  with  observations. 

From  iu^results'W-^concludei  that  the  non-uniform  cooling  of  the  post- 
flare corona  can  be  understood  as  a direct  consequence  of  the  temperature  and 
density  dependence  of  the  radiative  losses  from  a hlgh-temperature  solar  plasma. 
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